Actual source code: nepbasic.c

  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:    Note:
 42:    The default `NEP` type is `NEPRII`.

 44:    Level: beginner

 46: .seealso: [](ch:nep), `NEPSetUp()`, `NEPSolve()`, `NEPDestroy()`, `NEP`
 47: @*/
 48: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep)
 49: {
 50:   NEP            nep;

 52:   PetscFunctionBegin;
 53:   PetscAssertPointer(outnep,2);
 54:   PetscCall(NEPInitializePackage());
 55:   PetscCall(SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView));

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

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

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

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

123:   PetscCall(PetscNew(&nep->sc));
124:   *outnep = nep;
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: /*@
129:    NEPSetType - Selects the particular solver to be used in the `NEP` object.

131:    Logically Collective

133:    Input Parameters:
134: +  nep      - the nonlinear eigensolver context
135: -  type     - a known method

137:    Options Database Key:
138: .  -nep_type <method> - Sets the method; use `-help` for a list
139:     of available methods

141:    Notes:
142:    See `NEPType` for available methods. The default is `NEPRII`.

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

152:    Level: intermediate

154: .seealso: [](ch:nep), `NEPType`
155: @*/
156: PetscErrorCode NEPSetType(NEP nep,NEPType type)
157: {
158:   PetscErrorCode (*r)(NEP);
159:   PetscBool      match;

161:   PetscFunctionBegin;
163:   PetscAssertPointer(type,2);

165:   PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
166:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

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

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

174:   nep->state = NEP_STATE_INITIAL;
175:   PetscCall(PetscObjectChangeTypeName((PetscObject)nep,type));
176:   PetscCall((*r)(nep));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*@
181:    NEPGetType - Gets the `NEP` type as a string from the `NEP` object.

183:    Not Collective

185:    Input Parameter:
186: .  nep - the nonlinear eigensolver context

188:    Output Parameter:
189: .  type - name of `NEP` method

191:    Level: intermediate

193: .seealso: [](ch:nep), `NEPSetType()`
194: @*/
195: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
196: {
197:   PetscFunctionBegin;
199:   PetscAssertPointer(type,2);
200:   *type = ((PetscObject)nep)->type_name;
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

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

207:    Not Collective

209:    Input Parameters:
210: +  name - name of a new user-defined solver
211: -  function - routine to create the solver context

213:    Note:
214:    `NEPRegister()` may be called multiple times to add several user-defined solvers.

216:    Example Usage:
217: .vb
218:    NEPRegister("my_solver",MySolverCreate);
219: .ve

221:    Then, your solver can be chosen with the procedural interface via
222: .vb
223:    NEPSetType(nep,"my_solver")
224: .ve
225:    or at runtime via the option `-nep_type my_solver`.

227:    Level: advanced

229: .seealso: [](ch:nep), `NEPRegisterAll()`
230: @*/
231: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
232: {
233:   PetscFunctionBegin;
234:   PetscCall(NEPInitializePackage());
235:   PetscCall(PetscFunctionListAdd(&NEPList,name,function));
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: /*@C
240:    NEPMonitorRegister - Registers a `NEP` monitor routine that may be accessed with
241:    `NEPMonitorSetFromOptions()`.

243:    Not Collective

245:    Input Parameters:
246: +  name    - name of a new monitor routine
247: .  vtype   - a `PetscViewerType` for the output
248: .  format  - a `PetscViewerFormat` for the output
249: .  monitor - monitor routine, see `NEPMonitorRegisterFn`
250: .  create  - creation routine, or `NULL`
251: -  destroy - destruction routine, or `NULL`

253:    Notes:
254:    `NEPMonitorRegister()` may be called multiple times to add several user-defined monitors.

256:    The calling sequence for the given function matches the calling sequence of `NEPMonitorFn`
257:    functions passed to `NEPMonitorSet()` with the additional requirement that its final argument
258:    be a `PetscViewerAndFormat`.

260:    Example Usage:
261: .vb
262:    NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
263: .ve

265:    Then, your monitor can be chosen with the procedural interface via
266: .vb
267:    NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL);
268: .ve
269:    or at runtime via the option `-nep_monitor_my_monitor`.

271:    Level: advanced

273: .seealso: [](ch:nep), `NEPMonitorSet()`, `NEPMonitorRegisterAll()`, `NEPMonitorSetFromOptions()`
274: @*/
275: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
276: {
277:   char           key[PETSC_MAX_PATH_LEN];

279:   PetscFunctionBegin;
280:   PetscCall(NEPInitializePackage());
281:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
282:   PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
283:   if (create)  PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
284:   if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /*
289:    NEPReset_Problem - Destroys the problem matrices.
290: */
291: PetscErrorCode NEPReset_Problem(NEP nep)
292: {
293:   PetscInt       i;

295:   PetscFunctionBegin;
297:   PetscCall(MatDestroy(&nep->function));
298:   PetscCall(MatDestroy(&nep->function_pre));
299:   PetscCall(MatDestroy(&nep->jacobian));
300:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
301:     PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
302:     for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
303:     PetscCall(PetscFree(nep->f));
304:     PetscCall(PetscFree(nep->nrma));
305:     if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
306:     nep->nt = 0;
307:   }
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }
310: /*@
311:    NEPReset - Resets the NEP context to the initial state (prior to setup)
312:    and destroys any allocated Vecs and Mats.

314:    Collective

316:    Input Parameter:
317: .  nep - the nonlinear eigensolver context

319:    Level: advanced

321: .seealso: [](ch:nep), `NEPDestroy()`
322: @*/
323: PetscErrorCode NEPReset(NEP nep)
324: {
325:   PetscFunctionBegin;
327:   if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
328:   PetscTryTypeMethod(nep,reset);
329:   if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
330:   PetscCall(NEPReset_Problem(nep));
331:   PetscCall(BVDestroy(&nep->V));
332:   PetscCall(BVDestroy(&nep->W));
333:   PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
334:   PetscCall(MatDestroy(&nep->resolvent));
335:   nep->nwork = 0;
336:   nep->state = NEP_STATE_INITIAL;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: /*@
341:    NEPDestroy - Destroys the NEP context.

343:    Collective

345:    Input Parameter:
346: .  nep - the nonlinear eigensolver context

348:    Level: beginner

350: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPSolve()`
351: @*/
352: PetscErrorCode NEPDestroy(NEP *nep)
353: {
354:   PetscFunctionBegin;
355:   if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
357:   if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
358:   PetscCall(NEPReset(*nep));
359:   PetscTryTypeMethod(*nep,destroy);
360:   if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
361:   PetscCall(RGDestroy(&(*nep)->rg));
362:   PetscCall(DSDestroy(&(*nep)->ds));
363:   PetscCall(KSPDestroy(&(*nep)->refineksp));
364:   PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
365:   PetscCall(PetscFree((*nep)->sc));
366:   /* just in case the initial vectors have not been used */
367:   PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
368:   if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
369:   if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
370:   PetscCall(NEPMonitorCancel(*nep));
371:   PetscCall(PetscHeaderDestroy(nep));
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

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

378:    Collective

380:    Input Parameters:
381: +  nep - the nonlinear eigensolver context
382: -  bv  - the basis vectors object

384:    Note:
385:    Use NEPGetBV() to retrieve the basis vectors context (for example,
386:    to free it at the end of the computations).

388:    Level: advanced

390: .seealso: [](ch:nep), `NEPGetBV()`
391: @*/
392: PetscErrorCode NEPSetBV(NEP nep,BV bv)
393: {
394:   PetscFunctionBegin;
397:   PetscCheckSameComm(nep,1,bv,2);
398:   PetscCall(PetscObjectReference((PetscObject)bv));
399:   PetscCall(BVDestroy(&nep->V));
400:   nep->V = bv;
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: /*@
405:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
406:    eigensolver object.

408:    Not Collective

410:    Input Parameter:
411: .  nep - the nonlinear eigensolver context

413:    Output Parameter:
414: .  bv - basis vectors context

416:    Level: advanced

418: .seealso: [](ch:nep), `NEPSetBV()`
419: @*/
420: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
421: {
422:   PetscFunctionBegin;
424:   PetscAssertPointer(bv,2);
425:   if (!nep->V) {
426:     PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
427:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
428:     PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
429:   }
430:   *bv = nep->V;
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: /*@
435:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

437:    Collective

439:    Input Parameters:
440: +  nep - the nonlinear eigensolver context
441: -  rg  - the region object

443:    Note:
444:    Use NEPGetRG() to retrieve the region context (for example,
445:    to free it at the end of the computations).

447:    Level: advanced

449: .seealso: [](ch:nep), `NEPGetRG()`
450: @*/
451: PetscErrorCode NEPSetRG(NEP nep,RG rg)
452: {
453:   PetscFunctionBegin;
455:   if (rg) {
457:     PetscCheckSameComm(nep,1,rg,2);
458:   }
459:   PetscCall(PetscObjectReference((PetscObject)rg));
460:   PetscCall(RGDestroy(&nep->rg));
461:   nep->rg = rg;
462:   PetscFunctionReturn(PETSC_SUCCESS);
463: }

465: /*@
466:    NEPGetRG - Obtain the region object associated to the
467:    nonlinear eigensolver object.

469:    Not Collective

471:    Input Parameter:
472: .  nep - the nonlinear eigensolver context

474:    Output Parameter:
475: .  rg - region context

477:    Level: advanced

479: .seealso: [](ch:nep), `NEPSetRG()`
480: @*/
481: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
482: {
483:   PetscFunctionBegin;
485:   PetscAssertPointer(rg,2);
486:   if (!nep->rg) {
487:     PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
488:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
489:     PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
490:   }
491:   *rg = nep->rg;
492:   PetscFunctionReturn(PETSC_SUCCESS);
493: }

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

498:    Collective

500:    Input Parameters:
501: +  nep - the nonlinear eigensolver context
502: -  ds  - the direct solver object

504:    Note:
505:    Use NEPGetDS() to retrieve the direct solver context (for example,
506:    to free it at the end of the computations).

508:    Level: advanced

510: .seealso: [](ch:nep), `NEPGetDS()`
511: @*/
512: PetscErrorCode NEPSetDS(NEP nep,DS ds)
513: {
514:   PetscFunctionBegin;
517:   PetscCheckSameComm(nep,1,ds,2);
518:   PetscCall(PetscObjectReference((PetscObject)ds));
519:   PetscCall(DSDestroy(&nep->ds));
520:   nep->ds = ds;
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: /*@
525:    NEPGetDS - Obtain the direct solver object associated to the
526:    nonlinear eigensolver object.

528:    Not Collective

530:    Input Parameter:
531: .  nep - the nonlinear eigensolver context

533:    Output Parameter:
534: .  ds - direct solver context

536:    Level: advanced

538: .seealso: [](ch:nep), `NEPSetDS()`
539: @*/
540: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
541: {
542:   PetscFunctionBegin;
544:   PetscAssertPointer(ds,2);
545:   if (!nep->ds) {
546:     PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
547:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
548:     PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
549:   }
550:   *ds = nep->ds;
551:   PetscFunctionReturn(PETSC_SUCCESS);
552: }

554: /*@
555:    NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
556:    object in the refinement phase.

558:    Collective

560:    Input Parameter:
561: .  nep - the nonlinear eigensolver context

563:    Output Parameter:
564: .  ksp - ksp context

566:    Level: advanced

568: .seealso: [](ch:nep), `NEPSetRefine()`
569: @*/
570: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
571: {
572:   MPI_Comm       comm;

574:   PetscFunctionBegin;
576:   PetscAssertPointer(ksp,2);
577:   if (!nep->refineksp) {
578:     if (nep->npart>1) {
579:       /* Split in subcomunicators */
580:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
581:       PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
582:       PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
583:       PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
584:     } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
585:     PetscCall(KSPCreate(comm,&nep->refineksp));
586:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
587:     PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
588:     PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
589:     PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
590:     PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
591:   }
592:   *ksp = nep->refineksp;
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: /*@
597:    NEPSetTarget - Sets the value of the target.

599:    Logically Collective

601:    Input Parameters:
602: +  nep    - the nonlinear eigensolver context
603: -  target - the value of the target

605:    Options Database Key:
606: .  -nep_target <scalar> - the value of the target

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

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

616:    Level: intermediate

618: .seealso: [](ch:nep), `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
619: @*/
620: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
621: {
622:   PetscFunctionBegin;
625:   nep->target = target;
626:   PetscFunctionReturn(PETSC_SUCCESS);
627: }

629: /*@
630:    NEPGetTarget - Gets the value of the target.

632:    Not Collective

634:    Input Parameter:
635: .  nep - the nonlinear eigensolver context

637:    Output Parameter:
638: .  target - the value of the target

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

643:    Level: intermediate

645: .seealso: [](ch:nep), `NEPSetTarget()`
646: @*/
647: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
648: {
649:   PetscFunctionBegin;
651:   PetscAssertPointer(target,2);
652:   *target = nep->target;
653:   PetscFunctionReturn(PETSC_SUCCESS);
654: }

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

660:    Collective

662:    Input Parameters:
663: +  nep - the nonlinear eigensolver context
664: .  F   - Function matrix
665: .  P   - preconditioner matrix (usually same as F)
666: .  fun - Function evaluation routine (if NULL then NEP retains any
667:          previously set value), see NEPFunctionFn for the calling sequence
668: -  ctx - [optional] user-defined context for private data for the Function
669:          evaluation routine (may be NULL) (if NULL then NEP retains any
670:          previously set value)

672:    Level: beginner

674: .seealso: [](ch:nep), `NEPGetFunction()`, `NEPSetJacobian()`
675: @*/
676: PetscErrorCode NEPSetFunction(NEP nep,Mat F,Mat P,NEPFunctionFn *fun,void *ctx)
677: {
678:   PetscFunctionBegin;
682:   if (F) PetscCheckSameComm(nep,1,F,2);
683:   if (P) PetscCheckSameComm(nep,1,P,3);

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

688:   if (fun) nep->computefunction = fun;
689:   if (ctx) nep->functionctx     = ctx;
690:   if (F) {
691:     PetscCall(PetscObjectReference((PetscObject)F));
692:     PetscCall(MatDestroy(&nep->function));
693:     nep->function = F;
694:   }
695:   if (P) {
696:     PetscCall(PetscObjectReference((PetscObject)P));
697:     PetscCall(MatDestroy(&nep->function_pre));
698:     nep->function_pre = P;
699:   }
700:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
701:   nep->state = NEP_STATE_INITIAL;
702:   PetscFunctionReturn(PETSC_SUCCESS);
703: }

705: /*@C
706:    NEPGetFunction - Returns the Function matrix and optionally the user
707:    provided context for evaluating the Function.

709:    Collective

711:    Input Parameter:
712: .  nep - the nonlinear eigensolver context

714:    Output Parameters:
715: +  F   - location to stash Function matrix (or NULL)
716: .  P   - location to stash preconditioner matrix (or NULL)
717: .  fun - location to put Function function (or NULL)
718: -  ctx - location to stash Function context (or NULL)

720:    Level: advanced

722: .seealso: [](ch:nep), `NEPSetFunction()`
723: @*/
724: PetscErrorCode NEPGetFunction(NEP nep,Mat *F,Mat *P,NEPFunctionFn **fun,void **ctx)
725: {
726:   PetscFunctionBegin;
728:   NEPCheckCallback(nep,1);
729:   if (F)   *F   = nep->function;
730:   if (P)   *P   = nep->function_pre;
731:   if (fun) *fun = nep->computefunction;
732:   if (ctx) *ctx = nep->functionctx;
733:   PetscFunctionReturn(PETSC_SUCCESS);
734: }

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

740:    Collective

742:    Input Parameters:
743: +  nep - the nonlinear eigensolver context
744: .  J   - Jacobian matrix
745: .  jac - Jacobian evaluation routine (if NULL then NEP retains any
746:          previously set value), see NEPJacobianFn for the calling sequence
747: -  ctx - [optional] user-defined context for private data for the Jacobian
748:          evaluation routine (may be NULL) (if NULL then NEP retains any
749:          previously set value)

751:    Level: beginner

753: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetJacobian()`
754: @*/
755: PetscErrorCode NEPSetJacobian(NEP nep,Mat J,NEPJacobianFn *jac,void *ctx)
756: {
757:   PetscFunctionBegin;
760:   if (J) PetscCheckSameComm(nep,1,J,2);

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

765:   if (jac) nep->computejacobian = jac;
766:   if (ctx) nep->jacobianctx     = ctx;
767:   if (J) {
768:     PetscCall(PetscObjectReference((PetscObject)J));
769:     PetscCall(MatDestroy(&nep->jacobian));
770:     nep->jacobian = J;
771:   }
772:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
773:   nep->state = NEP_STATE_INITIAL;
774:   PetscFunctionReturn(PETSC_SUCCESS);
775: }

777: /*@C
778:    NEPGetJacobian - Returns the Jacobian matrix and optionally the user
779:    provided routine and context for evaluating the Jacobian.

781:    Collective

783:    Input Parameter:
784: .  nep - the nonlinear eigensolver context

786:    Output Parameters:
787: +  J   - location to stash Jacobian matrix (or NULL)
788: .  jac - location to put Jacobian function (or NULL)
789: -  ctx - location to stash Jacobian context (or NULL)

791:    Level: advanced

793: .seealso: [](ch:nep), `NEPSetJacobian()`
794: @*/
795: PetscErrorCode NEPGetJacobian(NEP nep,Mat *J,NEPJacobianFn **jac,void **ctx)
796: {
797:   PetscFunctionBegin;
799:   NEPCheckCallback(nep,1);
800:   if (J)   *J   = nep->jacobian;
801:   if (jac) *jac = nep->computejacobian;
802:   if (ctx) *ctx = nep->jacobianctx;
803:   PetscFunctionReturn(PETSC_SUCCESS);
804: }

806: /*@
807:    NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
808:    in split form.

810:    Collective

812:    Input Parameters:
813: +  nep - the nonlinear eigensolver context
814: .  nt  - number of terms in the split form
815: .  A   - array of matrices
816: .  f   - array of functions
817: -  str - structure flag for matrices

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

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

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

835:    Level: beginner

837: .seealso: [](ch:nep), `NEPGetSplitOperatorTerm()`, `NEPGetSplitOperatorInfo()`, `NEPSetSplitPreconditioner()`
838: @*/
839: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
840: {
841:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;

843:   PetscFunctionBegin;
846:   PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
847:   PetscAssertPointer(A,3);
848:   PetscAssertPointer(f,4);

851:   for (i=0;i<nt;i++) {
853:     PetscCheckSameComm(nep,1,A[i],3);
855:     PetscCheckSameComm(nep,1,f[i],4);
856:     PetscCall(MatGetSize(A[i],&m,&n));
857:     PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
858:     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);
859:     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);
860:     if (!i) { m0 = m; mloc0 = mloc; }
861:     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);
862:     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);
863:     PetscCall(PetscObjectReference((PetscObject)A[i]));
864:     PetscCall(PetscObjectReference((PetscObject)f[i]));
865:   }

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

870:   /* allocate space and copy matrices and functions */
871:   PetscCall(PetscMalloc1(nt,&nep->A));
872:   for (i=0;i<nt;i++) nep->A[i] = A[i];
873:   PetscCall(PetscMalloc1(nt,&nep->f));
874:   for (i=0;i<nt;i++) nep->f[i] = f[i];
875:   PetscCall(PetscCalloc1(nt,&nep->nrma));
876:   nep->nt    = nt;
877:   nep->mstr  = str;
878:   nep->fui   = NEP_USER_INTERFACE_SPLIT;
879:   nep->state = NEP_STATE_INITIAL;
880:   PetscFunctionReturn(PETSC_SUCCESS);
881: }

883: /*@
884:    NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
885:    the nonlinear operator in split form.

887:    Collective

889:    Input Parameters:
890: +  nep - the nonlinear eigensolver context
891: -  k   - the index of the requested term (starting in 0)

893:    Output Parameters:
894: +  A - the matrix of the requested term
895: -  f - the function of the requested term

897:    Level: intermediate

899: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorInfo()`
900: @*/
901: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
902: {
903:   PetscFunctionBegin;
906:   NEPCheckSplit(nep,1);
907:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
908:   if (A) *A = nep->A[k];
909:   if (f) *f = nep->f[k];
910:   PetscFunctionReturn(PETSC_SUCCESS);
911: }

913: /*@
914:    NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
915:    the nonlinear operator, as well as the structure flag for matrices.

917:    Not Collective

919:    Input Parameter:
920: .  nep - the nonlinear eigensolver context

922:    Output Parameters:
923: +  n   - the number of terms passed in NEPSetSplitOperator()
924: -  str - the matrix structure flag passed in NEPSetSplitOperator()

926:    Level: intermediate

928: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorTerm()`
929: @*/
930: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
931: {
932:   PetscFunctionBegin;
934:   NEPCheckSplit(nep,1);
935:   if (n)   *n = nep->nt;
936:   if (str) *str = nep->mstr;
937:   PetscFunctionReturn(PETSC_SUCCESS);
938: }

940: /*@
941:    NEPSetSplitPreconditioner - Sets an operator in split form from which
942:    to build the preconditioner to be used when solving the nonlinear
943:    eigenvalue problem in split form.

945:    Collective

947:    Input Parameters:
948: +  nep  - the nonlinear eigensolver context
949: .  ntp  - number of terms in the split preconditioner
950: .  P    - array of matrices
951: -  strp - structure flag for matrices

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

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

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

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

970:    Level: advanced

972: .seealso: [](ch:nep), `NEPGetSplitPreconditionerTerm()`, `NEPGetSplitPreconditionerInfo()`, `NEPSetSplitOperator()`
973: @*/
974: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
975: {
976:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;

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

987:   for (i=0;i<ntp;i++) {
989:     PetscCheckSameComm(nep,1,P[i],3);
990:     PetscCall(MatGetSize(P[i],&m,&n));
991:     PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
992:     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);
993:     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);
994:     if (!i) { m0 = m; mloc0 = mloc; }
995:     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);
996:     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);
997:     PetscCall(PetscObjectReference((PetscObject)P[i]));
998:   }

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

1003:   /* allocate space and copy matrices */
1004:   if (ntp) {
1005:     PetscCall(PetscMalloc1(ntp,&nep->P));
1006:     for (i=0;i<ntp;i++) nep->P[i] = P[i];
1007:   }
1008:   nep->mstrp = strp;
1009:   nep->state = NEP_STATE_INITIAL;
1010:   PetscFunctionReturn(PETSC_SUCCESS);
1011: }

1013: /*@
1014:    NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1015:    the split preconditioner.

1017:    Not Collective

1019:    Input Parameters:
1020: +  nep - the nonlinear eigensolver context
1021: -  k   - the index of the requested term (starting in 0)

1023:    Output Parameter:
1024: .  P  - the matrix of the requested term

1026:    Level: advanced

1028: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerInfo()`
1029: @*/
1030: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1031: {
1032:   PetscFunctionBegin;
1035:   PetscAssertPointer(P,3);
1036:   NEPCheckSplit(nep,1);
1037:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1038:   PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1039:   *P = nep->P[k];
1040:   PetscFunctionReturn(PETSC_SUCCESS);
1041: }

1043: /*@
1044:    NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1045:    preconditioner, as well as the structure flag for matrices.

1047:    Not Collective

1049:    Input Parameter:
1050: .  nep - the nonlinear eigensolver context

1052:    Output Parameters:
1053: +  n    - the number of terms passed in NEPSetSplitPreconditioner()
1054: -  strp - the matrix structure flag passed in NEPSetSplitPreconditioner()

1056:    Level: advanced

1058: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerTerm()`
1059: @*/
1060: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1061: {
1062:   PetscFunctionBegin;
1064:   NEPCheckSplit(nep,1);
1065:   if (n)    *n    = nep->P? nep->nt: 0;
1066:   if (strp) *strp = nep->mstrp;
1067:   PetscFunctionReturn(PETSC_SUCCESS);
1068: }