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: Level: intermediate
192: .seealso: [](ch:nep), `NEPSetType()`
193: @*/
194: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
195: {
196: PetscFunctionBegin;
198: PetscAssertPointer(type,2);
199: *type = ((PetscObject)nep)->type_name;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@C
204: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
206: Not Collective
208: Input Parameters:
209: + name - name of a new user-defined solver
210: - function - routine to create the solver context
212: Note:
213: `NEPRegister()` may be called multiple times to add several user-defined solvers.
215: Example Usage:
216: .vb
217: NEPRegister("my_solver",MySolverCreate);
218: .ve
220: Then, your solver can be chosen with the procedural interface via
221: .vb
222: NEPSetType(nep,"my_solver")
223: .ve
224: or at runtime via the option `-nep_type my_solver`.
226: Level: advanced
228: .seealso: [](ch:nep), `NEPRegisterAll()`
229: @*/
230: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
231: {
232: PetscFunctionBegin;
233: PetscCall(NEPInitializePackage());
234: PetscCall(PetscFunctionListAdd(&NEPList,name,function));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@C
239: NEPMonitorRegister - Registers a `NEP` monitor routine that may be accessed with
240: `NEPMonitorSetFromOptions()`.
242: Not Collective
244: Input Parameters:
245: + name - name of a new monitor routine
246: . vtype - a `PetscViewerType` for the output
247: . format - a `PetscViewerFormat` for the output
248: . monitor - monitor routine, see `NEPMonitorRegisterFn`
249: . create - creation routine, or `NULL`
250: - destroy - destruction routine, or `NULL`
252: Notes:
253: `NEPMonitorRegister()` may be called multiple times to add several user-defined monitors.
255: The calling sequence for the given function matches the calling sequence of `NEPMonitorFn`
256: functions passed to `NEPMonitorSet()` with the additional requirement that its final argument
257: be a `PetscViewerAndFormat`.
259: Example Usage:
260: .vb
261: NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
262: .ve
264: Then, your monitor can be chosen with the procedural interface via
265: .vb
266: NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL);
267: .ve
268: or at runtime via the option `-nep_monitor_my_monitor`.
270: Level: advanced
272: .seealso: [](ch:nep), `NEPMonitorSet()`, `NEPMonitorRegisterAll()`, `NEPMonitorSetFromOptions()`
273: @*/
274: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
275: {
276: char key[PETSC_MAX_PATH_LEN];
278: PetscFunctionBegin;
279: PetscCall(NEPInitializePackage());
280: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
281: PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
282: if (create) PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
283: if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*
288: NEPReset_Problem - Destroys the problem matrices.
289: */
290: PetscErrorCode NEPReset_Problem(NEP nep)
291: {
292: PetscInt i;
294: PetscFunctionBegin;
296: PetscCall(MatDestroy(&nep->function));
297: PetscCall(MatDestroy(&nep->function_pre));
298: PetscCall(MatDestroy(&nep->jacobian));
299: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
300: PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
301: for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
302: PetscCall(PetscFree(nep->f));
303: PetscCall(PetscFree(nep->nrma));
304: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
305: nep->nt = 0;
306: }
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
309: /*@
310: NEPReset - Resets the `NEP` context to the initial state (prior to setup)
311: and destroys any allocated `Vec`s and `Mat`s.
313: Collective
315: Input Parameter:
316: . nep - the nonlinear eigensolver context
318: Level: advanced
320: .seealso: [](ch:nep), `NEPDestroy()`
321: @*/
322: PetscErrorCode NEPReset(NEP nep)
323: {
324: PetscFunctionBegin;
326: if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
327: PetscTryTypeMethod(nep,reset);
328: if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
329: PetscCall(NEPReset_Problem(nep));
330: PetscCall(BVDestroy(&nep->V));
331: PetscCall(BVDestroy(&nep->W));
332: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
333: PetscCall(MatDestroy(&nep->resolvent));
334: nep->nwork = 0;
335: nep->state = NEP_STATE_INITIAL;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*@
340: NEPDestroy - Destroys the `NEP` context.
342: Collective
344: Input Parameter:
345: . nep - the nonlinear eigensolver context
347: Level: beginner
349: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPSolve()`
350: @*/
351: PetscErrorCode NEPDestroy(NEP *nep)
352: {
353: PetscFunctionBegin;
354: if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
356: if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
357: PetscCall(NEPReset(*nep));
358: PetscTryTypeMethod(*nep,destroy);
359: if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
360: PetscCall(RGDestroy(&(*nep)->rg));
361: PetscCall(DSDestroy(&(*nep)->ds));
362: PetscCall(KSPDestroy(&(*nep)->refineksp));
363: PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
364: PetscCall(PetscFree((*nep)->sc));
365: /* just in case the initial vectors have not been used */
366: PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
367: if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
368: if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
369: PetscCall(NEPMonitorCancel(*nep));
370: PetscCall(PetscHeaderDestroy(nep));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*@
375: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
377: Collective
379: Input Parameters:
380: + nep - the nonlinear eigensolver context
381: - bv - the basis vectors object
383: Note:
384: Use `NEPGetBV()` to retrieve the basis vectors context (for example,
385: to free it at the end of the computations).
387: Level: advanced
389: .seealso: [](ch:nep), `NEPGetBV()`
390: @*/
391: PetscErrorCode NEPSetBV(NEP nep,BV bv)
392: {
393: PetscFunctionBegin;
396: PetscCheckSameComm(nep,1,bv,2);
397: PetscCall(PetscObjectReference((PetscObject)bv));
398: PetscCall(BVDestroy(&nep->V));
399: nep->V = bv;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@
404: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
405: eigensolver object.
407: Not Collective
409: Input Parameter:
410: . nep - the nonlinear eigensolver context
412: Output Parameter:
413: . bv - basis vectors context
415: Level: advanced
417: .seealso: [](ch:nep), `NEPSetBV()`
418: @*/
419: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
420: {
421: PetscFunctionBegin;
423: PetscAssertPointer(bv,2);
424: if (!nep->V) {
425: PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
426: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
427: PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
428: }
429: *bv = nep->V;
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@
434: NEPSetRG - Associates a region object to the nonlinear eigensolver.
436: Collective
438: Input Parameters:
439: + nep - the nonlinear eigensolver context
440: - rg - the region object
442: Note:
443: Use `NEPGetRG()` to retrieve the region context (for example,
444: to free it at the end of the computations).
446: Level: advanced
448: .seealso: [](ch:nep), `NEPGetRG()`
449: @*/
450: PetscErrorCode NEPSetRG(NEP nep,RG rg)
451: {
452: PetscFunctionBegin;
454: if (rg) {
456: PetscCheckSameComm(nep,1,rg,2);
457: }
458: PetscCall(PetscObjectReference((PetscObject)rg));
459: PetscCall(RGDestroy(&nep->rg));
460: nep->rg = rg;
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: NEPGetRG - Obtain the region object associated to the
466: nonlinear eigensolver object.
468: Not Collective
470: Input Parameter:
471: . nep - the nonlinear eigensolver context
473: Output Parameter:
474: . rg - region context
476: Level: advanced
478: .seealso: [](ch:nep), `NEPSetRG()`
479: @*/
480: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
481: {
482: PetscFunctionBegin;
484: PetscAssertPointer(rg,2);
485: if (!nep->rg) {
486: PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
487: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
488: PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
489: }
490: *rg = nep->rg;
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
497: Collective
499: Input Parameters:
500: + nep - the nonlinear eigensolver context
501: - ds - the direct solver object
503: Note:
504: Use `NEPGetDS()` to retrieve the direct solver context (for example,
505: to free it at the end of the computations).
507: Level: advanced
509: .seealso: [](ch:nep), `NEPGetDS()`
510: @*/
511: PetscErrorCode NEPSetDS(NEP nep,DS ds)
512: {
513: PetscFunctionBegin;
516: PetscCheckSameComm(nep,1,ds,2);
517: PetscCall(PetscObjectReference((PetscObject)ds));
518: PetscCall(DSDestroy(&nep->ds));
519: nep->ds = ds;
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: NEPGetDS - Obtain the direct solver object associated to the
525: nonlinear eigensolver object.
527: Not Collective
529: Input Parameter:
530: . nep - the nonlinear eigensolver context
532: Output Parameter:
533: . ds - direct solver context
535: Level: advanced
537: .seealso: [](ch:nep), `NEPSetDS()`
538: @*/
539: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
540: {
541: PetscFunctionBegin;
543: PetscAssertPointer(ds,2);
544: if (!nep->ds) {
545: PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
546: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
547: PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
548: }
549: *ds = nep->ds;
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: /*@
554: NEPRefineGetKSP - Obtain the `KSP` object used by the eigensolver
555: object in the refinement phase.
557: Collective
559: Input Parameter:
560: . nep - the nonlinear eigensolver context
562: Output Parameter:
563: . ksp - the linear solver context
565: Level: advanced
567: .seealso: [](ch:nep), `NEPSetRefine()`
568: @*/
569: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
570: {
571: MPI_Comm comm;
573: PetscFunctionBegin;
575: PetscAssertPointer(ksp,2);
576: if (!nep->refineksp) {
577: if (nep->npart>1) {
578: /* Split in subcomunicators */
579: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
580: PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
581: PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
582: PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
583: } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
584: PetscCall(KSPCreate(comm,&nep->refineksp));
585: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
586: PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
587: PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
588: PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
589: PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
590: }
591: *ksp = nep->refineksp;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: /*@
596: NEPSetTarget - Sets the value of the target.
598: Logically Collective
600: Input Parameters:
601: + nep - the nonlinear eigensolver context
602: - target - the value of the target
604: Options Database Key:
605: . -nep_target \<target\> - the value of the target
607: Notes:
608: The target is a scalar value used to determine the portion of the spectrum
609: of interest. It is used in combination with `NEPSetWhichEigenpairs()`.
611: When PETSc is built with real scalars, it is not possible to specify a
612: complex target.
614: In the case of complex scalars, a complex value can be provided in the
615: command line with `[+/-][realnumber][+/-]realnumberi` with no spaces, e.g.
616: `-nep_target 1.0+2.0i`.
618: Level: intermediate
620: .seealso: [](ch:nep), `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
621: @*/
622: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
623: {
624: PetscFunctionBegin;
627: nep->target = target;
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*@
632: NEPGetTarget - Gets the value of the target.
634: Not Collective
636: Input Parameter:
637: . nep - the nonlinear eigensolver context
639: Output Parameter:
640: . target - the value of the target
642: Note:
643: If the target was not set by the user, then zero is returned.
645: Level: intermediate
647: .seealso: [](ch:nep), `NEPSetTarget()`
648: @*/
649: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
650: {
651: PetscFunctionBegin;
653: PetscAssertPointer(target,2);
654: *target = nep->target;
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*@C
659: NEPSetFunction - Sets the function to compute the nonlinear Function $T(\lambda)$
660: as well as the location to store the matrix.
662: Collective
664: Input Parameters:
665: + nep - the nonlinear eigensolver context
666: . F - Function matrix
667: . P - preconditioner matrix (usually the same as `F`)
668: . fun - Function evaluation routine (if `NULL` then `NEP` retains any
669: previously set value), see `NEPFunctionFn` for the calling sequence
670: - ctx - [optional] user-defined context for private data for the Function
671: evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
672: previously set value)
674: Level: beginner
676: .seealso: [](ch:nep), `NEPGetFunction()`, `NEPSetJacobian()`
677: @*/
678: PetscErrorCode NEPSetFunction(NEP nep,Mat F,Mat P,NEPFunctionFn *fun,void *ctx)
679: {
680: PetscFunctionBegin;
684: if (F) PetscCheckSameComm(nep,1,F,2);
685: if (P) PetscCheckSameComm(nep,1,P,3);
687: if (nep->state) PetscCall(NEPReset(nep));
688: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
690: if (fun) nep->computefunction = fun;
691: if (ctx) nep->functionctx = ctx;
692: if (F) {
693: PetscCall(PetscObjectReference((PetscObject)F));
694: PetscCall(MatDestroy(&nep->function));
695: nep->function = F;
696: }
697: if (P) {
698: PetscCall(PetscObjectReference((PetscObject)P));
699: PetscCall(MatDestroy(&nep->function_pre));
700: nep->function_pre = P;
701: }
702: nep->fui = NEP_USER_INTERFACE_CALLBACK;
703: nep->state = NEP_STATE_INITIAL;
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@C
708: NEPGetFunction - Returns the Function matrix and optionally the user
709: provided context for evaluating the Function.
711: Collective
713: Input Parameter:
714: . nep - the nonlinear eigensolver context
716: Output Parameters:
717: + F - location to stash Function matrix (or `NULL`)
718: . P - location to stash preconditioner matrix (or `NULL`)
719: . fun - location to put Function function (or `NULL`)
720: - ctx - location to stash Function context (or `NULL`)
722: Level: advanced
724: .seealso: [](ch:nep), `NEPSetFunction()`
725: @*/
726: PetscErrorCode NEPGetFunction(NEP nep,Mat *F,Mat *P,NEPFunctionFn **fun,void **ctx)
727: {
728: PetscFunctionBegin;
730: NEPCheckCallback(nep,1);
731: if (F) *F = nep->function;
732: if (P) *P = nep->function_pre;
733: if (fun) *fun = nep->computefunction;
734: if (ctx) *ctx = nep->functionctx;
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: /*@C
739: NEPSetJacobian - Sets the function to compute the Jacobian $T'(\lambda)$ as well
740: as the location to store the matrix.
742: Collective
744: Input Parameters:
745: + nep - the nonlinear eigensolver context
746: . J - Jacobian matrix
747: . jac - Jacobian evaluation routine (if `NULL` then `NEP` retains any
748: previously set value), see `NEPJacobianFn` for the calling sequence
749: - ctx - [optional] user-defined context for private data for the Jacobian
750: evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
751: previously set value)
753: Level: beginner
755: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetJacobian()`
756: @*/
757: PetscErrorCode NEPSetJacobian(NEP nep,Mat J,NEPJacobianFn *jac,void *ctx)
758: {
759: PetscFunctionBegin;
762: if (J) PetscCheckSameComm(nep,1,J,2);
764: if (nep->state) PetscCall(NEPReset(nep));
765: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
767: if (jac) nep->computejacobian = jac;
768: if (ctx) nep->jacobianctx = ctx;
769: if (J) {
770: PetscCall(PetscObjectReference((PetscObject)J));
771: PetscCall(MatDestroy(&nep->jacobian));
772: nep->jacobian = J;
773: }
774: nep->fui = NEP_USER_INTERFACE_CALLBACK;
775: nep->state = NEP_STATE_INITIAL;
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: /*@C
780: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
781: provided routine and context for evaluating the Jacobian.
783: Collective
785: Input Parameter:
786: . nep - the nonlinear eigensolver context
788: Output Parameters:
789: + J - location to stash Jacobian matrix (or `NULL`)
790: . jac - location to put Jacobian function (or `NULL`)
791: - ctx - location to stash Jacobian context (or `NULL`)
793: Level: advanced
795: .seealso: [](ch:nep), `NEPSetJacobian()`
796: @*/
797: PetscErrorCode NEPGetJacobian(NEP nep,Mat *J,NEPJacobianFn **jac,void **ctx)
798: {
799: PetscFunctionBegin;
801: NEPCheckCallback(nep,1);
802: if (J) *J = nep->jacobian;
803: if (jac) *jac = nep->computejacobian;
804: if (ctx) *ctx = nep->jacobianctx;
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: /*@
809: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
810: in split form.
812: Collective
814: Input Parameters:
815: + nep - the nonlinear eigensolver context
816: . nt - number of terms in the split form
817: . A - array of matrices
818: . f - array of functions
819: - str - structure flag for matrices
821: Notes:
822: The nonlinear operator is written as $T(\lambda) = \sum_i A_i f_i(\lambda)$,
823: for $i=1,\dots,n$. The derivative $T'(\lambda)$ can be obtained using the
824: derivatives of $f_i$.
826: The structure flag provides information about $A_i$'s nonzero pattern
827: (see `MatStructure`). If all matrices have the same pattern, then
828: use `SAME_NONZERO_PATTERN`. If the patterns are different but contained
829: in the pattern of the first one, then use `SUBSET_NONZERO_PATTERN`. If
830: patterns are known to be different, use `DIFFERENT_NONZERO_PATTERN`.
831: If set to `UNKNOWN_NONZERO_PATTERN`, the patterns will be compared to
832: determine if they are equal.
834: This function must be called before `NEPSetUp()`. If it is called again
835: after `NEPSetUp()` then the `NEP` object is reset.
837: Level: beginner
839: .seealso: [](ch:nep), `NEPGetSplitOperatorTerm()`, `NEPGetSplitOperatorInfo()`, `NEPSetSplitPreconditioner()`
840: @*/
841: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
842: {
843: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
845: PetscFunctionBegin;
848: PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
849: PetscAssertPointer(A,3);
850: PetscAssertPointer(f,4);
853: for (i=0;i<nt;i++) {
855: PetscCheckSameComm(nep,1,A[i],3);
857: PetscCheckSameComm(nep,1,f[i],4);
858: PetscCall(MatGetSize(A[i],&m,&n));
859: PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
860: 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);
861: 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);
862: if (!i) { m0 = m; mloc0 = mloc; }
863: 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);
864: 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);
865: PetscCall(PetscObjectReference((PetscObject)A[i]));
866: PetscCall(PetscObjectReference((PetscObject)f[i]));
867: }
869: if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
870: else PetscCall(NEPReset_Problem(nep));
872: /* allocate space and copy matrices and functions */
873: PetscCall(PetscMalloc1(nt,&nep->A));
874: for (i=0;i<nt;i++) nep->A[i] = A[i];
875: PetscCall(PetscMalloc1(nt,&nep->f));
876: for (i=0;i<nt;i++) nep->f[i] = f[i];
877: PetscCall(PetscCalloc1(nt,&nep->nrma));
878: nep->nt = nt;
879: nep->mstr = str;
880: nep->fui = NEP_USER_INTERFACE_SPLIT;
881: nep->state = NEP_STATE_INITIAL;
882: PetscFunctionReturn(PETSC_SUCCESS);
883: }
885: /*@
886: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
887: the nonlinear operator in split form.
889: Collective
891: Input Parameters:
892: + nep - the nonlinear eigensolver context
893: - k - the index of the requested term (starting in 0)
895: Output Parameters:
896: + A - the matrix of the requested term
897: - f - the function of the requested term
899: Level: intermediate
901: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorInfo()`
902: @*/
903: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
904: {
905: PetscFunctionBegin;
908: NEPCheckSplit(nep,1);
909: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
910: if (A) *A = nep->A[k];
911: if (f) *f = nep->f[k];
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*@
916: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
917: the nonlinear operator, as well as the structure flag for matrices.
919: Not Collective
921: Input Parameter:
922: . nep - the nonlinear eigensolver context
924: Output Parameters:
925: + n - the number of terms passed in `NEPSetSplitOperator()`
926: - str - the matrix structure flag passed in `NEPSetSplitOperator()`
928: Level: intermediate
930: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorTerm()`
931: @*/
932: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
933: {
934: PetscFunctionBegin;
936: NEPCheckSplit(nep,1);
937: if (n) *n = nep->nt;
938: if (str) *str = nep->mstr;
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: /*@
943: NEPSetSplitPreconditioner - Sets an operator in split form from which
944: to build the preconditioner to be used when solving the nonlinear
945: eigenvalue problem in split form.
947: Collective
949: Input Parameters:
950: + nep - the nonlinear eigensolver context
951: . ntp - number of terms in the split preconditioner
952: . P - array of matrices
953: - strp - structure flag for matrices
955: Notes:
956: The matrix for the preconditioner is expressed as $P(\lambda) =
957: \sum_i P_i f_i(\lambda)$, for $i=1,\dots,n$, where the $f_i$ functions
958: are the same as in `NEPSetSplitOperator()`. It is not necessary to call
959: this function. If it is not invoked, then the preconditioner is
960: built from $T(\lambda)$, i.e., both matrices and functions passed in
961: `NEPSetSplitOperator()`.
963: The structure flag provides information about $P_i$'s nonzero pattern
964: in the same way as in `NEPSetSplitOperator()`.
966: If the functions defining the preconditioner operator were different
967: from the ones given in `NEPSetSplitOperator()`, then the split form
968: cannot be used. Use the callback interface instead.
970: Use `ntp=0` to reset a previously set split preconditioner.
972: Level: advanced
974: .seealso: [](ch:nep), `NEPGetSplitPreconditionerTerm()`, `NEPGetSplitPreconditionerInfo()`, `NEPSetSplitOperator()`
975: @*/
976: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
977: {
978: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
980: PetscFunctionBegin;
983: PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
984: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
985: PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
986: if (ntp) PetscAssertPointer(P,3);
989: for (i=0;i<ntp;i++) {
991: PetscCheckSameComm(nep,1,P[i],3);
992: PetscCall(MatGetSize(P[i],&m,&n));
993: PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
994: 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);
995: 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);
996: if (!i) { m0 = m; mloc0 = mloc; }
997: 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);
998: 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);
999: PetscCall(PetscObjectReference((PetscObject)P[i]));
1000: }
1002: PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
1003: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
1005: /* allocate space and copy matrices */
1006: if (ntp) {
1007: PetscCall(PetscMalloc1(ntp,&nep->P));
1008: for (i=0;i<ntp;i++) nep->P[i] = P[i];
1009: }
1010: nep->mstrp = strp;
1011: nep->state = NEP_STATE_INITIAL;
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: /*@
1016: NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1017: the split preconditioner.
1019: Not Collective
1021: Input Parameters:
1022: + nep - the nonlinear eigensolver context
1023: - k - the index of the requested term (starting in 0)
1025: Output Parameter:
1026: . P - the matrix of the requested term
1028: Level: advanced
1030: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerInfo()`
1031: @*/
1032: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1033: {
1034: PetscFunctionBegin;
1037: PetscAssertPointer(P,3);
1038: NEPCheckSplit(nep,1);
1039: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1040: PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1041: *P = nep->P[k];
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@
1046: NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1047: preconditioner, as well as the structure flag for matrices.
1049: Not Collective
1051: Input Parameter:
1052: . nep - the nonlinear eigensolver context
1054: Output Parameters:
1055: + n - the number of terms passed in `NEPSetSplitPreconditioner()`
1056: - strp - the matrix structure flag passed in `NEPSetSplitPreconditioner()`
1058: Level: advanced
1060: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerTerm()`
1061: @*/
1062: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1063: {
1064: PetscFunctionBegin;
1066: NEPCheckSplit(nep,1);
1067: if (n) *n = nep->P? nep->nt: 0;
1068: if (strp) *strp = nep->mstrp;
1069: PetscFunctionReturn(PETSC_SUCCESS);
1070: }