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 `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 type - sets the method; use `-help` for a list of available methods
140: Notes:
141: See `NEPType` for available methods. The default is `NEPRII`.
143: Normally, it is best to use the `NEPSetFromOptions()` command and
144: then set the `NEP` type from the options database rather than by using
145: this routine. Using the options database provides the user with
146: maximum flexibility in evaluating the different available methods.
147: The `NEPSetType()` routine is provided for those situations where it
148: is necessary to set the iterative solver independently of the command
149: line or options database.
151: Level: intermediate
153: .seealso: [](ch:nep), `NEPType`
154: @*/
155: PetscErrorCode NEPSetType(NEP nep,NEPType type)
156: {
157: PetscErrorCode (*r)(NEP);
158: PetscBool match;
160: PetscFunctionBegin;
162: PetscAssertPointer(type,2);
164: PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
165: if (match) PetscFunctionReturn(PETSC_SUCCESS);
167: PetscCall(PetscFunctionListFind(NEPList,type,&r));
168: PetscCheck(r,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);
170: PetscTryTypeMethod(nep,destroy);
171: PetscCall(PetscMemzero(nep->ops,sizeof(struct _NEPOps)));
173: nep->state = NEP_STATE_INITIAL;
174: PetscCall(PetscObjectChangeTypeName((PetscObject)nep,type));
175: PetscCall((*r)(nep));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: NEPGetType - Gets the `NEP` type as a string from the `NEP` object.
182: Not Collective
184: Input Parameter:
185: . nep - the nonlinear eigensolver context
187: Output Parameter:
188: . type - name of `NEP` method
190: Note:
191: `type` should not be retained for later use as it will be an invalid pointer
192: if the `NEPType` of `nep` is changed.
194: Level: intermediate
196: .seealso: [](ch:nep), `NEPSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
197: @*/
198: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
199: {
200: PetscFunctionBegin;
202: PetscAssertPointer(type,2);
203: *type = ((PetscObject)nep)->type_name;
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@C
208: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
210: Not Collective
212: Input Parameters:
213: + name - name of a new user-defined solver
214: - function - routine to create the solver context
216: Note:
217: `NEPRegister()` may be called multiple times to add several user-defined solvers.
219: Example Usage:
220: .vb
221: NEPRegister("my_solver",MySolverCreate);
222: .ve
224: Then, your solver can be chosen with the procedural interface via
225: .vb
226: NEPSetType(nep,"my_solver")
227: .ve
228: or at runtime via the option `-nep_type my_solver`.
230: Level: advanced
232: .seealso: [](ch:nep), `NEPRegisterAll()`
233: @*/
234: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
235: {
236: PetscFunctionBegin;
237: PetscCall(NEPInitializePackage());
238: PetscCall(PetscFunctionListAdd(&NEPList,name,function));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@C
243: NEPMonitorRegister - Registers a `NEP` monitor routine that may be accessed with
244: `NEPMonitorSetFromOptions()`.
246: Not Collective
248: Input Parameters:
249: + name - name of a new monitor routine
250: . vtype - a `PetscViewerType` for the output
251: . format - a `PetscViewerFormat` for the output
252: . monitor - monitor routine, see `NEPMonitorRegisterFn`
253: . create - creation routine, or `NULL`
254: - destroy - destruction routine, or `NULL`
256: Notes:
257: `NEPMonitorRegister()` may be called multiple times to add several user-defined monitors.
259: The calling sequence for the given function matches the calling sequence of `NEPMonitorFn`
260: functions passed to `NEPMonitorSet()` with the additional requirement that its final argument
261: be a `PetscViewerAndFormat`.
263: Example Usage:
264: .vb
265: NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
266: .ve
268: Then, your monitor can be chosen with the procedural interface via
269: .vb
270: NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL);
271: .ve
272: or at runtime via the option `-nep_monitor_my_monitor`.
274: Level: advanced
276: .seealso: [](ch:nep), `NEPMonitorSet()`, `NEPMonitorRegisterAll()`, `NEPMonitorSetFromOptions()`
277: @*/
278: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
279: {
280: char key[PETSC_MAX_PATH_LEN];
282: PetscFunctionBegin;
283: PetscCall(NEPInitializePackage());
284: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
285: PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
286: if (create) PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
287: if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*
292: NEPReset_Problem - Destroys the problem matrices.
293: */
294: PetscErrorCode NEPReset_Problem(NEP nep)
295: {
296: PetscInt i;
298: PetscFunctionBegin;
300: PetscCall(MatDestroy(&nep->function));
301: PetscCall(MatDestroy(&nep->function_pre));
302: PetscCall(MatDestroy(&nep->jacobian));
303: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
304: PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
305: for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
306: PetscCall(PetscFree(nep->f));
307: PetscCall(PetscFree(nep->nrma));
308: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
309: nep->nt = 0;
310: }
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
313: /*@
314: NEPReset - Resets the `NEP` context to the initial state (prior to setup)
315: and destroys any allocated `Vec`s and `Mat`s.
317: Collective
319: Input Parameter:
320: . nep - the nonlinear eigensolver context
322: Level: advanced
324: .seealso: [](ch:nep), `NEPDestroy()`
325: @*/
326: PetscErrorCode NEPReset(NEP nep)
327: {
328: PetscFunctionBegin;
330: if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
331: PetscTryTypeMethod(nep,reset);
332: if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
333: PetscCall(NEPReset_Problem(nep));
334: PetscCall(BVDestroy(&nep->V));
335: PetscCall(BVDestroy(&nep->W));
336: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
337: PetscCall(MatDestroy(&nep->resolvent));
338: nep->nwork = 0;
339: nep->state = NEP_STATE_INITIAL;
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@
344: NEPDestroy - Destroys the `NEP` context.
346: Collective
348: Input Parameter:
349: . nep - the nonlinear eigensolver context
351: Level: beginner
353: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPSolve()`
354: @*/
355: PetscErrorCode NEPDestroy(NEP *nep)
356: {
357: PetscFunctionBegin;
358: if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
360: if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
361: PetscCall(NEPReset(*nep));
362: PetscTryTypeMethod(*nep,destroy);
363: if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
364: PetscCall(RGDestroy(&(*nep)->rg));
365: PetscCall(DSDestroy(&(*nep)->ds));
366: PetscCall(KSPDestroy(&(*nep)->refineksp));
367: PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
368: PetscCall(PetscFree((*nep)->sc));
369: /* just in case the initial vectors have not been used */
370: PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
371: if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
372: if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
373: PetscCall(NEPMonitorCancel(*nep));
374: PetscCall(PetscHeaderDestroy(nep));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
381: Collective
383: Input Parameters:
384: + nep - the nonlinear eigensolver context
385: - bv - the basis vectors object
387: Note:
388: Use `NEPGetBV()` to retrieve the basis vectors context (for example,
389: to free it at the end of the computations).
391: Level: advanced
393: .seealso: [](ch:nep), `NEPGetBV()`
394: @*/
395: PetscErrorCode NEPSetBV(NEP nep,BV bv)
396: {
397: PetscFunctionBegin;
400: PetscCheckSameComm(nep,1,bv,2);
401: PetscCall(PetscObjectReference((PetscObject)bv));
402: PetscCall(BVDestroy(&nep->V));
403: nep->V = bv;
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: /*@
408: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
409: eigensolver object.
411: Not Collective
413: Input Parameter:
414: . nep - the nonlinear eigensolver context
416: Output Parameter:
417: . bv - basis vectors context
419: Level: advanced
421: .seealso: [](ch:nep), `NEPSetBV()`
422: @*/
423: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
424: {
425: PetscFunctionBegin;
427: PetscAssertPointer(bv,2);
428: if (!nep->V) {
429: PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
430: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
431: PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
432: }
433: *bv = nep->V;
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: /*@
438: NEPSetRG - Associates a region object to the nonlinear eigensolver.
440: Collective
442: Input Parameters:
443: + nep - the nonlinear eigensolver context
444: - rg - the region object
446: Note:
447: Use `NEPGetRG()` to retrieve the region context (for example,
448: to free it at the end of the computations).
450: Level: advanced
452: .seealso: [](ch:nep), `NEPGetRG()`
453: @*/
454: PetscErrorCode NEPSetRG(NEP nep,RG rg)
455: {
456: PetscFunctionBegin;
458: if (rg) {
460: PetscCheckSameComm(nep,1,rg,2);
461: }
462: PetscCall(PetscObjectReference((PetscObject)rg));
463: PetscCall(RGDestroy(&nep->rg));
464: nep->rg = rg;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: NEPGetRG - Obtain the region object associated to the
470: nonlinear eigensolver object.
472: Not Collective
474: Input Parameter:
475: . nep - the nonlinear eigensolver context
477: Output Parameter:
478: . rg - region context
480: Level: advanced
482: .seealso: [](ch:nep), `NEPSetRG()`
483: @*/
484: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
485: {
486: PetscFunctionBegin;
488: PetscAssertPointer(rg,2);
489: if (!nep->rg) {
490: PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
491: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
492: PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
493: }
494: *rg = nep->rg;
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
501: Collective
503: Input Parameters:
504: + nep - the nonlinear eigensolver context
505: - ds - the direct solver object
507: Note:
508: Use `NEPGetDS()` to retrieve the direct solver context (for example,
509: to free it at the end of the computations).
511: Level: advanced
513: .seealso: [](ch:nep), `NEPGetDS()`
514: @*/
515: PetscErrorCode NEPSetDS(NEP nep,DS ds)
516: {
517: PetscFunctionBegin;
520: PetscCheckSameComm(nep,1,ds,2);
521: PetscCall(PetscObjectReference((PetscObject)ds));
522: PetscCall(DSDestroy(&nep->ds));
523: nep->ds = ds;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: NEPGetDS - Obtain the direct solver object associated to the
529: nonlinear eigensolver object.
531: Not Collective
533: Input Parameter:
534: . nep - the nonlinear eigensolver context
536: Output Parameter:
537: . ds - direct solver context
539: Level: advanced
541: .seealso: [](ch:nep), `NEPSetDS()`
542: @*/
543: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
544: {
545: PetscFunctionBegin;
547: PetscAssertPointer(ds,2);
548: if (!nep->ds) {
549: PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
550: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
551: PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
552: }
553: *ds = nep->ds;
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: /*@
558: NEPRefineGetKSP - Obtain the `KSP` object used by the eigensolver
559: object in the refinement phase.
561: Collective
563: Input Parameter:
564: . nep - the nonlinear eigensolver context
566: Output Parameter:
567: . ksp - the linear solver context
569: Level: advanced
571: .seealso: [](ch:nep), `NEPSetRefine()`
572: @*/
573: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
574: {
575: MPI_Comm comm;
577: PetscFunctionBegin;
579: PetscAssertPointer(ksp,2);
580: if (!nep->refineksp) {
581: if (nep->npart>1) {
582: /* Split in subcomunicators */
583: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
584: PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
585: PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
586: PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
587: } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
588: PetscCall(KSPCreate(comm,&nep->refineksp));
589: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
590: PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
591: PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
592: PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
593: PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
594: }
595: *ksp = nep->refineksp;
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*@
600: NEPSetTarget - Sets the value of the target.
602: Logically Collective
604: Input Parameters:
605: + nep - the nonlinear eigensolver context
606: - target - the value of the target
608: Options Database Key:
609: . -nep_target target - the value of the target
611: Notes:
612: The target is a scalar value used to determine the portion of the spectrum
613: of interest. It is used in combination with `NEPSetWhichEigenpairs()`.
615: When PETSc is built with real scalars, it is not possible to specify a
616: complex target.
618: In the case of complex scalars, a complex value can be provided in the
619: command line with `[+/-][realnumber][+/-]realnumberi` with no spaces, e.g.
620: `-nep_target 1.0+2.0i`.
622: Level: intermediate
624: .seealso: [](ch:nep), `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
625: @*/
626: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
627: {
628: PetscFunctionBegin;
631: nep->target = target;
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*@
636: NEPGetTarget - Gets the value of the target.
638: Not Collective
640: Input Parameter:
641: . nep - the nonlinear eigensolver context
643: Output Parameter:
644: . target - the value of the target
646: Note:
647: If the target was not set by the user, then zero is returned.
649: Level: intermediate
651: .seealso: [](ch:nep), `NEPSetTarget()`
652: @*/
653: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
654: {
655: PetscFunctionBegin;
657: PetscAssertPointer(target,2);
658: *target = nep->target;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@C
663: NEPSetFunction - Sets the function to compute the nonlinear Function $T(\lambda)$
664: as well as the location to store the matrix.
666: Collective
668: Input Parameters:
669: + nep - the nonlinear eigensolver context
670: . F - Function matrix
671: . P - preconditioner matrix (usually the same as `F`)
672: . fun - Function evaluation routine (if `NULL` then `NEP` retains any
673: previously set value), see `NEPFunctionFn` for the calling sequence
674: - ctx - [optional] user-defined context for private data for the Function
675: evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
676: previously set value)
678: Level: beginner
680: .seealso: [](ch:nep), `NEPGetFunction()`, `NEPSetJacobian()`
681: @*/
682: PetscErrorCode NEPSetFunction(NEP nep,Mat F,Mat P,NEPFunctionFn *fun,PetscCtx ctx)
683: {
684: PetscFunctionBegin;
688: if (F) PetscCheckSameComm(nep,1,F,2);
689: if (P) PetscCheckSameComm(nep,1,P,3);
691: if (nep->state) PetscCall(NEPReset(nep));
692: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
694: if (fun) nep->computefunction = fun;
695: if (ctx) nep->functionctx = ctx;
696: if (F) {
697: PetscCall(PetscObjectReference((PetscObject)F));
698: PetscCall(MatDestroy(&nep->function));
699: nep->function = F;
700: }
701: if (P) {
702: PetscCall(PetscObjectReference((PetscObject)P));
703: PetscCall(MatDestroy(&nep->function_pre));
704: nep->function_pre = P;
705: }
706: nep->fui = NEP_USER_INTERFACE_CALLBACK;
707: nep->state = NEP_STATE_INITIAL;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@C
712: NEPGetFunction - Returns the Function matrix and optionally the user
713: provided context for evaluating the Function.
715: Collective
717: Input Parameter:
718: . nep - the nonlinear eigensolver context
720: Output Parameters:
721: + F - location to stash Function matrix (or `NULL`)
722: . P - location to stash preconditioner matrix (or `NULL`)
723: . fun - location to put Function function (or `NULL`)
724: - ctx - location to stash Function context (or `NULL`)
726: Level: advanced
728: .seealso: [](ch:nep), `NEPSetFunction()`
729: @*/
730: PetscErrorCode NEPGetFunction(NEP nep,Mat *F,Mat *P,NEPFunctionFn **fun,PetscCtxRt ctx)
731: {
732: PetscFunctionBegin;
734: NEPCheckCallback(nep,1);
735: if (F) *F = nep->function;
736: if (P) *P = nep->function_pre;
737: if (fun) *fun = nep->computefunction;
738: if (ctx) *(void**)ctx = nep->functionctx;
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: /*@C
743: NEPSetJacobian - Sets the function to compute the Jacobian $T'(\lambda)$ as well
744: as the location to store the matrix.
746: Collective
748: Input Parameters:
749: + nep - the nonlinear eigensolver context
750: . J - Jacobian matrix
751: . jac - Jacobian evaluation routine (if `NULL` then `NEP` retains any
752: previously set value), see `NEPJacobianFn` for the calling sequence
753: - ctx - [optional] user-defined context for private data for the Jacobian
754: evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
755: previously set value)
757: Level: beginner
759: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetJacobian()`
760: @*/
761: PetscErrorCode NEPSetJacobian(NEP nep,Mat J,NEPJacobianFn *jac,PetscCtx ctx)
762: {
763: PetscFunctionBegin;
766: if (J) PetscCheckSameComm(nep,1,J,2);
768: if (nep->state) PetscCall(NEPReset(nep));
769: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
771: if (jac) nep->computejacobian = jac;
772: if (ctx) nep->jacobianctx = ctx;
773: if (J) {
774: PetscCall(PetscObjectReference((PetscObject)J));
775: PetscCall(MatDestroy(&nep->jacobian));
776: nep->jacobian = J;
777: }
778: nep->fui = NEP_USER_INTERFACE_CALLBACK;
779: nep->state = NEP_STATE_INITIAL;
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: /*@C
784: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
785: provided routine and context for evaluating the Jacobian.
787: Collective
789: Input Parameter:
790: . nep - the nonlinear eigensolver context
792: Output Parameters:
793: + J - location to stash Jacobian matrix (or `NULL`)
794: . jac - location to put Jacobian function (or `NULL`)
795: - ctx - location to stash Jacobian context (or `NULL`)
797: Level: advanced
799: .seealso: [](ch:nep), `NEPSetJacobian()`
800: @*/
801: PetscErrorCode NEPGetJacobian(NEP nep,Mat *J,NEPJacobianFn **jac,PetscCtxRt ctx)
802: {
803: PetscFunctionBegin;
805: NEPCheckCallback(nep,1);
806: if (J) *J = nep->jacobian;
807: if (jac) *jac = nep->computejacobian;
808: if (ctx) *(void**)ctx = nep->jacobianctx;
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: /*@
813: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
814: in split form.
816: Collective
818: Input Parameters:
819: + nep - the nonlinear eigensolver context
820: . nt - number of terms in the split form
821: . A - array of matrices
822: . f - array of functions
823: - str - structure flag for matrices
825: Notes:
826: The nonlinear operator is written as $T(\lambda) = \sum_i A_i f_i(\lambda)$,
827: for $i=1,\dots,n$. The derivative $T'(\lambda)$ can be obtained using the
828: derivatives of $f_i$.
830: The structure flag provides information about $A_i$'s nonzero pattern
831: (see `MatStructure`). If all matrices have the same pattern, then
832: use `SAME_NONZERO_PATTERN`. If the patterns are different but contained
833: in the pattern of the first one, then use `SUBSET_NONZERO_PATTERN`. If
834: patterns are known to be different, use `DIFFERENT_NONZERO_PATTERN`.
835: If set to `UNKNOWN_NONZERO_PATTERN`, the patterns will be compared to
836: determine if they are equal.
838: This function must be called before `NEPSetUp()`. If it is called again
839: after `NEPSetUp()` then the `NEP` object is reset.
841: Level: beginner
843: .seealso: [](ch:nep), `NEPGetSplitOperatorTerm()`, `NEPGetSplitOperatorInfo()`, `NEPSetSplitPreconditioner()`
844: @*/
845: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
846: {
847: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
849: PetscFunctionBegin;
852: PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
853: PetscAssertPointer(A,3);
854: PetscAssertPointer(f,4);
857: for (i=0;i<nt;i++) {
859: PetscCheckSameComm(nep,1,A[i],3);
861: PetscCheckSameComm(nep,1,f[i],4);
862: PetscCall(MatGetSize(A[i],&m,&n));
863: PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
864: 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);
865: 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);
866: if (!i) { m0 = m; mloc0 = mloc; }
867: 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);
868: 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);
869: PetscCall(PetscObjectReference((PetscObject)A[i]));
870: PetscCall(PetscObjectReference((PetscObject)f[i]));
871: }
873: if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
874: else PetscCall(NEPReset_Problem(nep));
876: /* allocate space and copy matrices and functions */
877: PetscCall(PetscMalloc1(nt,&nep->A));
878: for (i=0;i<nt;i++) nep->A[i] = A[i];
879: PetscCall(PetscMalloc1(nt,&nep->f));
880: for (i=0;i<nt;i++) nep->f[i] = f[i];
881: PetscCall(PetscCalloc1(nt,&nep->nrma));
882: nep->nt = nt;
883: nep->mstr = str;
884: nep->fui = NEP_USER_INTERFACE_SPLIT;
885: nep->state = NEP_STATE_INITIAL;
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: /*@
890: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
891: the nonlinear operator in split form.
893: Collective
895: Input Parameters:
896: + nep - the nonlinear eigensolver context
897: - k - the index of the requested term (starting in 0)
899: Output Parameters:
900: + A - the matrix of the requested term
901: - f - the function of the requested term
903: Level: intermediate
905: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorInfo()`
906: @*/
907: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
908: {
909: PetscFunctionBegin;
912: NEPCheckSplit(nep,1);
913: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
914: if (A) *A = nep->A[k];
915: if (f) *f = nep->f[k];
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /*@
920: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
921: the nonlinear operator, as well as the structure flag for matrices.
923: Not Collective
925: Input Parameter:
926: . nep - the nonlinear eigensolver context
928: Output Parameters:
929: + n - the number of terms passed in `NEPSetSplitOperator()`
930: - str - the matrix structure flag passed in `NEPSetSplitOperator()`
932: Level: intermediate
934: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorTerm()`
935: @*/
936: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
937: {
938: PetscFunctionBegin;
940: NEPCheckSplit(nep,1);
941: if (n) *n = nep->nt;
942: if (str) *str = nep->mstr;
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: NEPSetSplitPreconditioner - Sets an operator in split form from which
948: to build the preconditioner to be used when solving the nonlinear
949: eigenvalue problem in split form.
951: Collective
953: Input Parameters:
954: + nep - the nonlinear eigensolver context
955: . ntp - number of terms in the split preconditioner
956: . P - array of matrices
957: - strp - structure flag for matrices
959: Notes:
960: The matrix for the preconditioner is expressed as $P(\lambda) =
961: \sum_i P_i f_i(\lambda)$, for $i=1,\dots,n$, where the $f_i$ functions
962: are the same as in `NEPSetSplitOperator()`. It is not necessary to call
963: this function. If it is not invoked, then the preconditioner is
964: built from $T(\lambda)$, i.e., both matrices and functions passed in
965: `NEPSetSplitOperator()`.
967: The structure flag provides information about $P_i$'s nonzero pattern
968: in the same way as in `NEPSetSplitOperator()`.
970: If the functions defining the preconditioner operator were different
971: from the ones given in `NEPSetSplitOperator()`, then the split form
972: cannot be used. Use the callback interface instead.
974: Use `ntp=0` to reset a previously set split preconditioner.
976: Level: advanced
978: .seealso: [](ch:nep), `NEPGetSplitPreconditionerTerm()`, `NEPGetSplitPreconditionerInfo()`, `NEPSetSplitOperator()`
979: @*/
980: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
981: {
982: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
984: PetscFunctionBegin;
987: PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
988: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
989: PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
990: if (ntp) PetscAssertPointer(P,3);
993: for (i=0;i<ntp;i++) {
995: PetscCheckSameComm(nep,1,P[i],3);
996: PetscCall(MatGetSize(P[i],&m,&n));
997: PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
998: 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);
999: 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);
1000: if (!i) { m0 = m; mloc0 = mloc; }
1001: 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);
1002: 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);
1003: PetscCall(PetscObjectReference((PetscObject)P[i]));
1004: }
1006: PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
1007: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
1009: /* allocate space and copy matrices */
1010: if (ntp) {
1011: PetscCall(PetscMalloc1(ntp,&nep->P));
1012: for (i=0;i<ntp;i++) nep->P[i] = P[i];
1013: }
1014: nep->mstrp = strp;
1015: nep->state = NEP_STATE_INITIAL;
1016: PetscFunctionReturn(PETSC_SUCCESS);
1017: }
1019: /*@
1020: NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1021: the split preconditioner.
1023: Not Collective
1025: Input Parameters:
1026: + nep - the nonlinear eigensolver context
1027: - k - the index of the requested term (starting in 0)
1029: Output Parameter:
1030: . P - the matrix of the requested term
1032: Level: advanced
1034: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerInfo()`
1035: @*/
1036: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1037: {
1038: PetscFunctionBegin;
1041: PetscAssertPointer(P,3);
1042: NEPCheckSplit(nep,1);
1043: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1044: PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1045: *P = nep->P[k];
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: /*@
1050: NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1051: preconditioner, as well as the structure flag for matrices.
1053: Not Collective
1055: Input Parameter:
1056: . nep - the nonlinear eigensolver context
1058: Output Parameters:
1059: + n - the number of terms passed in `NEPSetSplitPreconditioner()`
1060: - strp - the matrix structure flag passed in `NEPSetSplitPreconditioner()`
1062: Level: advanced
1064: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerTerm()`
1065: @*/
1066: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1067: {
1068: PetscFunctionBegin;
1070: NEPCheckSplit(nep,1);
1071: if (n) *n = nep->P? nep->nt: 0;
1072: if (strp) *strp = nep->mstrp;
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }