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:    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 - Registers a  NEP monitor routine that may be accessed with NEPMonitorSetFromOptions().

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, see NEPMonitorRegisterFn
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:    The calling sequence for the given function matches the calling sequence of NEPMonitorFn
252:    functions passed to NEPMonitorSet() with the additional requirement that its final argument
253:    be a PetscViewerAndFormat.

255:    Example Usage:
256: .vb
257:    NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
258: .ve

260:    Then, your monitor can be chosen with the procedural interface via
261: $      NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL)
262:    or at runtime via the option
263: $      -nep_monitor_my_monitor

265:    Level: advanced

267: .seealso: `NEPMonitorSet()`, `NEPMonitorRegisterAll()`
268: @*/
269: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
270: {
271:   char           key[PETSC_MAX_PATH_LEN];

273:   PetscFunctionBegin;
274:   PetscCall(NEPInitializePackage());
275:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
276:   PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
277:   if (create)  PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
278:   if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: /*
283:    NEPReset_Problem - Destroys the problem matrices.
284: */
285: PetscErrorCode NEPReset_Problem(NEP nep)
286: {
287:   PetscInt       i;

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

308:    Collective

310:    Input Parameter:
311: .  nep - eigensolver context obtained from NEPCreate()

313:    Level: advanced

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

334: /*@
335:    NEPDestroy - Destroys the NEP context.

337:    Collective

339:    Input Parameter:
340: .  nep - eigensolver context obtained from NEPCreate()

342:    Level: beginner

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

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

372:    Collective

374:    Input Parameters:
375: +  nep - eigensolver context obtained from NEPCreate()
376: -  bv  - the basis vectors object

378:    Note:
379:    Use NEPGetBV() to retrieve the basis vectors context (for example,
380:    to free it at the end of the computations).

382:    Level: advanced

384: .seealso: `NEPGetBV()`
385: @*/
386: PetscErrorCode NEPSetBV(NEP nep,BV bv)
387: {
388:   PetscFunctionBegin;
391:   PetscCheckSameComm(nep,1,bv,2);
392:   PetscCall(PetscObjectReference((PetscObject)bv));
393:   PetscCall(BVDestroy(&nep->V));
394:   nep->V = bv;
395:   PetscFunctionReturn(PETSC_SUCCESS);
396: }

398: /*@
399:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
400:    eigensolver object.

402:    Not Collective

404:    Input Parameters:
405: .  nep - eigensolver context obtained from NEPCreate()

407:    Output Parameter:
408: .  bv - basis vectors context

410:    Level: advanced

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

428: /*@
429:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

431:    Collective

433:    Input Parameters:
434: +  nep - eigensolver context obtained from NEPCreate()
435: -  rg  - the region object

437:    Note:
438:    Use NEPGetRG() to retrieve the region context (for example,
439:    to free it at the end of the computations).

441:    Level: advanced

443: .seealso: `NEPGetRG()`
444: @*/
445: PetscErrorCode NEPSetRG(NEP nep,RG rg)
446: {
447:   PetscFunctionBegin;
449:   if (rg) {
451:     PetscCheckSameComm(nep,1,rg,2);
452:   }
453:   PetscCall(PetscObjectReference((PetscObject)rg));
454:   PetscCall(RGDestroy(&nep->rg));
455:   nep->rg = rg;
456:   PetscFunctionReturn(PETSC_SUCCESS);
457: }

459: /*@
460:    NEPGetRG - Obtain the region object associated to the
461:    nonlinear eigensolver object.

463:    Not Collective

465:    Input Parameters:
466: .  nep - eigensolver context obtained from NEPCreate()

468:    Output Parameter:
469: .  rg - region context

471:    Level: advanced

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

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

492:    Collective

494:    Input Parameters:
495: +  nep - eigensolver context obtained from NEPCreate()
496: -  ds  - the direct solver object

498:    Note:
499:    Use NEPGetDS() to retrieve the direct solver context (for example,
500:    to free it at the end of the computations).

502:    Level: advanced

504: .seealso: `NEPGetDS()`
505: @*/
506: PetscErrorCode NEPSetDS(NEP nep,DS ds)
507: {
508:   PetscFunctionBegin;
511:   PetscCheckSameComm(nep,1,ds,2);
512:   PetscCall(PetscObjectReference((PetscObject)ds));
513:   PetscCall(DSDestroy(&nep->ds));
514:   nep->ds = ds;
515:   PetscFunctionReturn(PETSC_SUCCESS);
516: }

518: /*@
519:    NEPGetDS - Obtain the direct solver object associated to the
520:    nonlinear eigensolver object.

522:    Not Collective

524:    Input Parameters:
525: .  nep - eigensolver context obtained from NEPCreate()

527:    Output Parameter:
528: .  ds - direct solver context

530:    Level: advanced

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

548: /*@
549:    NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
550:    object in the refinement phase.

552:    Collective

554:    Input Parameters:
555: .  nep - eigensolver context obtained from NEPCreate()

557:    Output Parameter:
558: .  ksp - ksp context

560:    Level: advanced

562: .seealso: `NEPSetRefine()`
563: @*/
564: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
565: {
566:   MPI_Comm       comm;

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

590: /*@
591:    NEPSetTarget - Sets the value of the target.

593:    Logically Collective

595:    Input Parameters:
596: +  nep    - eigensolver context
597: -  target - the value of the target

599:    Options Database Key:
600: .  -nep_target <scalar> - the value of the target

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

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

610:    Level: intermediate

612: .seealso: `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
613: @*/
614: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
615: {
616:   PetscFunctionBegin;
619:   nep->target = target;
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: /*@
624:    NEPGetTarget - Gets the value of the target.

626:    Not Collective

628:    Input Parameter:
629: .  nep - eigensolver context

631:    Output Parameter:
632: .  target - the value of the target

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

637:    Level: intermediate

639: .seealso: `NEPSetTarget()`
640: @*/
641: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
642: {
643:   PetscFunctionBegin;
645:   PetscAssertPointer(target,2);
646:   *target = nep->target;
647:   PetscFunctionReturn(PETSC_SUCCESS);
648: }

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

654:    Collective

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

666:    Level: beginner

668: .seealso: `NEPGetFunction()`, `NEPSetJacobian()`
669: @*/
670: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,NEPFunctionFn *fun,void *ctx)
671: {
672:   PetscFunctionBegin;
676:   if (A) PetscCheckSameComm(nep,1,A,2);
677:   if (B) PetscCheckSameComm(nep,1,B,3);

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

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

699: /*@C
700:    NEPGetFunction - Returns the Function matrix and optionally the user
701:    provided context for evaluating the Function.

703:    Collective

705:    Input Parameter:
706: .  nep - the nonlinear eigensolver context

708:    Output Parameters:
709: +  A   - location to stash Function matrix (or NULL)
710: .  B   - location to stash preconditioner matrix (or NULL)
711: .  fun - location to put Function function (or NULL)
712: -  ctx - location to stash Function context (or NULL)

714:    Level: advanced

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

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

734:    Collective

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

745:    Level: beginner

747: .seealso: `NEPSetFunction()`, `NEPGetJacobian()`
748: @*/
749: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,NEPJacobianFn *jac,void *ctx)
750: {
751:   PetscFunctionBegin;
754:   if (A) PetscCheckSameComm(nep,1,A,2);

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

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

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

775:    Collective

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

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

785:    Level: advanced

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

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

804:    Collective

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

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

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

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

829:    Level: beginner

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

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

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

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

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

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

881:    Collective

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

887:    Output Parameters:
888: +  A - the matrix of the requested term
889: -  f - the function of the requested term

891:    Level: intermediate

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

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

911:    Not Collective

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

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

920:    Level: intermediate

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

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

939:    Collective

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

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

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

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

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

964:    Level: advanced

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

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

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

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

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

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

1011:    Not Collective

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

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

1020:    Level: advanced

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

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

1041:    Not Collective

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

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

1050:    Level: advanced

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