Actual source code: nepbasic.c

slepc-3.22.1 2024-10-28
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:   PetscCall(NEPMonitorCancel(*nep));
360:   PetscCall(PetscHeaderDestroy(nep));
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

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

367:    Collective

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

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

377:    Level: advanced

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

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

397:    Not Collective

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

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

405:    Level: advanced

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

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

426:    Collective

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

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

436:    Level: advanced

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

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

458:    Not Collective

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

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

466:    Level: advanced

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

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

487:    Collective

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

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

497:    Level: advanced

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

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

517:    Not Collective

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

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

525:    Level: advanced

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

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

547:    Collective

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

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

555:    Level: advanced

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

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

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

588:    Logically Collective

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

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

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

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

605:    Level: intermediate

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

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

621:    Not Collective

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

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

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

632:    Level: intermediate

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

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

649:    Collective

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

661:    Level: beginner

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

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

677:   if (fun) nep->computefunction = fun;
678:   if (ctx) nep->functionctx     = ctx;
679:   if (A) {
680:     PetscCall(PetscObjectReference((PetscObject)A));
681:     PetscCall(MatDestroy(&nep->function));
682:     nep->function = A;
683:   }
684:   if (B) {
685:     PetscCall(PetscObjectReference((PetscObject)B));
686:     PetscCall(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(PETSC_SUCCESS);
692: }

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

698:    Not Collective

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,NEPFunctionFn **fun,void **ctx)
714: {
715:   PetscFunctionBegin;
717:   NEPCheckCallback(nep,1);
718:   if (A)   *A   = nep->function;
719:   if (B)   *B   = nep->function_pre;
720:   if (fun) *fun = nep->computefunction;
721:   if (ctx) *ctx = nep->functionctx;
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

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

729:    Collective

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

740:    Level: beginner

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

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

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

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

770:    Not Collective

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

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

780:    Level: advanced

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

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

799:    Collective

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

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

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

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

824:    Level: beginner

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

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

840:   for (i=0;i<nt;i++) {
842:     PetscCheckSameComm(nep,1,A[i],3);
844:     PetscCheckSameComm(nep,1,f[i],4);
845:     PetscCall(MatGetSize(A[i],&m,&n));
846:     PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
847:     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);
848:     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);
849:     if (!i) { m0 = m; mloc0 = mloc; }
850:     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);
851:     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);
852:     PetscCall(PetscObjectReference((PetscObject)A[i]));
853:     PetscCall(PetscObjectReference((PetscObject)f[i]));
854:   }

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

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

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

876:    Not Collective

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

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

886:    Level: intermediate

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

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

906:    Not Collective

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

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

915:    Level: intermediate

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

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

934:    Collective

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

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

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

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

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

959:    Level: advanced

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

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

976:   for (i=0;i<ntp;i++) {
978:     PetscCheckSameComm(nep,1,P[i],3);
979:     PetscCall(MatGetSize(P[i],&m,&n));
980:     PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
981:     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);
982:     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);
983:     if (!i) { m0 = m; mloc0 = mloc; }
984:     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);
985:     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);
986:     PetscCall(PetscObjectReference((PetscObject)P[i]));
987:   }

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

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

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

1006:    Not Collective

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

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

1015:    Level: advanced

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

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

1036:    Not Collective

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

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

1045:    Level: advanced

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