Actual source code: nepbasic.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Basic NEP routines
12: */
14: #include <slepc/private/nepimpl.h>
16: /* Logging support */
17: PetscClassId NEP_CLASSID = 0;
18: PetscLogEvent NEP_SetUp = 0,NEP_Solve = 0,NEP_Refine = 0,NEP_FunctionEval = 0,NEP_JacobianEval = 0,NEP_Resolvent = 0,NEP_CISS_SVD = 0;
20: /* List of registered NEP routines */
21: PetscFunctionList NEPList = NULL;
22: PetscBool NEPRegisterAllCalled = PETSC_FALSE;
24: /* List of registered NEP monitors */
25: PetscFunctionList NEPMonitorList = NULL;
26: PetscFunctionList NEPMonitorCreateList = NULL;
27: PetscFunctionList NEPMonitorDestroyList = NULL;
28: PetscBool NEPMonitorRegisterAllCalled = PETSC_FALSE;
30: /*@
31: NEPCreate - Creates the default `NEP` context.
33: Collective
35: Input Parameter:
36: . comm - MPI communicator
38: Output Parameter:
39: . outnep - location to put the `NEP` context
41: Note:
42: The default `NEP` type is `NEPRII`.
44: Level: beginner
46: .seealso: [](ch:nep), `NEPSetUp()`, `NEPSolve()`, `NEPDestroy()`, `NEP`
47: @*/
48: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep)
49: {
50: NEP nep;
52: PetscFunctionBegin;
53: PetscAssertPointer(outnep,2);
54: PetscCall(NEPInitializePackage());
55: PetscCall(SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView));
57: nep->max_it = PETSC_DETERMINE;
58: nep->nev = 1;
59: nep->ncv = PETSC_DETERMINE;
60: nep->mpd = PETSC_DETERMINE;
61: nep->nini = 0;
62: nep->target = 0.0;
63: nep->tol = PETSC_DETERMINE;
64: nep->conv = NEP_CONV_REL;
65: nep->stop = NEP_STOP_BASIC;
66: nep->which = (NEPWhich)0;
67: nep->problem_type = (NEPProblemType)0;
68: nep->refine = NEP_REFINE_NONE;
69: nep->npart = 1;
70: nep->rtol = PETSC_DETERMINE;
71: nep->rits = PETSC_DETERMINE;
72: nep->scheme = (NEPRefineScheme)0;
73: nep->trackall = PETSC_FALSE;
74: nep->twosided = PETSC_FALSE;
76: nep->computefunction = NULL;
77: nep->computejacobian = NULL;
78: nep->functionctx = NULL;
79: nep->jacobianctx = NULL;
80: nep->converged = NEPConvergedRelative;
81: nep->convergeduser = NULL;
82: nep->convergeddestroy= NULL;
83: nep->stopping = NEPStoppingBasic;
84: nep->stoppinguser = NULL;
85: nep->stoppingdestroy = NULL;
86: nep->convergedctx = NULL;
87: nep->stoppingctx = NULL;
88: nep->numbermonitors = 0;
90: nep->ds = NULL;
91: nep->V = NULL;
92: nep->W = NULL;
93: nep->rg = NULL;
94: nep->function = NULL;
95: nep->function_pre = NULL;
96: nep->jacobian = NULL;
97: nep->A = NULL;
98: nep->f = NULL;
99: nep->nt = 0;
100: nep->mstr = UNKNOWN_NONZERO_PATTERN;
101: nep->P = NULL;
102: nep->mstrp = UNKNOWN_NONZERO_PATTERN;
103: nep->IS = NULL;
104: nep->eigr = NULL;
105: nep->eigi = NULL;
106: nep->errest = NULL;
107: nep->perm = NULL;
108: nep->nwork = 0;
109: nep->work = NULL;
110: nep->data = NULL;
112: nep->state = NEP_STATE_INITIAL;
113: nep->nconv = 0;
114: nep->its = 0;
115: nep->n = 0;
116: nep->nloc = 0;
117: nep->nrma = NULL;
118: nep->fui = (NEPUserInterface)0;
119: nep->useds = PETSC_FALSE;
120: nep->resolvent = NULL;
121: nep->reason = NEP_CONVERGED_ITERATING;
123: PetscCall(PetscNew(&nep->sc));
124: *outnep = nep;
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*@
129: NEPSetType - Selects the particular solver to be used in the `NEP` object.
131: Logically Collective
133: Input Parameters:
134: + nep - the nonlinear eigensolver context
135: - type - a known method
137: Options Database Key:
138: . -nep_type <method> - Sets the method; use `-help` for a list
139: of available methods
141: Notes:
142: See `NEPType` for available methods. The default is `NEPRII`.
144: Normally, it is best to use the `NEPSetFromOptions()` command and
145: then set the `NEP` type from the options database rather than by using
146: this routine. Using the options database provides the user with
147: maximum flexibility in evaluating the different available methods.
148: The `NEPSetType()` routine is provided for those situations where it
149: is necessary to set the iterative solver independently of the command
150: line or options database.
152: Level: intermediate
154: .seealso: [](ch:nep), `NEPType`
155: @*/
156: PetscErrorCode NEPSetType(NEP nep,NEPType type)
157: {
158: PetscErrorCode (*r)(NEP);
159: PetscBool match;
161: PetscFunctionBegin;
163: PetscAssertPointer(type,2);
165: PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
166: if (match) PetscFunctionReturn(PETSC_SUCCESS);
168: PetscCall(PetscFunctionListFind(NEPList,type,&r));
169: PetscCheck(r,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);
171: PetscTryTypeMethod(nep,destroy);
172: PetscCall(PetscMemzero(nep->ops,sizeof(struct _NEPOps)));
174: nep->state = NEP_STATE_INITIAL;
175: PetscCall(PetscObjectChangeTypeName((PetscObject)nep,type));
176: PetscCall((*r)(nep));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: NEPGetType - Gets the `NEP` type as a string from the `NEP` object.
183: Not Collective
185: Input Parameter:
186: . nep - the nonlinear eigensolver context
188: Output Parameter:
189: . type - name of `NEP` method
191: Level: intermediate
193: .seealso: [](ch:nep), `NEPSetType()`
194: @*/
195: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
196: {
197: PetscFunctionBegin;
199: PetscAssertPointer(type,2);
200: *type = ((PetscObject)nep)->type_name;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@C
205: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
207: Not Collective
209: Input Parameters:
210: + name - name of a new user-defined solver
211: - function - routine to create the solver context
213: Note:
214: `NEPRegister()` may be called multiple times to add several user-defined solvers.
216: Example Usage:
217: .vb
218: NEPRegister("my_solver",MySolverCreate);
219: .ve
221: Then, your solver can be chosen with the procedural interface via
222: .vb
223: NEPSetType(nep,"my_solver")
224: .ve
225: or at runtime via the option `-nep_type my_solver`.
227: Level: advanced
229: .seealso: [](ch:nep), `NEPRegisterAll()`
230: @*/
231: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
232: {
233: PetscFunctionBegin;
234: PetscCall(NEPInitializePackage());
235: PetscCall(PetscFunctionListAdd(&NEPList,name,function));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@C
240: NEPMonitorRegister - Registers a `NEP` monitor routine that may be accessed with
241: `NEPMonitorSetFromOptions()`.
243: Not Collective
245: Input Parameters:
246: + name - name of a new monitor routine
247: . vtype - a `PetscViewerType` for the output
248: . format - a `PetscViewerFormat` for the output
249: . monitor - monitor routine, see `NEPMonitorRegisterFn`
250: . create - creation routine, or `NULL`
251: - destroy - destruction routine, or `NULL`
253: Notes:
254: `NEPMonitorRegister()` may be called multiple times to add several user-defined monitors.
256: The calling sequence for the given function matches the calling sequence of `NEPMonitorFn`
257: functions passed to `NEPMonitorSet()` with the additional requirement that its final argument
258: be a `PetscViewerAndFormat`.
260: Example Usage:
261: .vb
262: NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
263: .ve
265: Then, your monitor can be chosen with the procedural interface via
266: .vb
267: NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL);
268: .ve
269: or at runtime via the option `-nep_monitor_my_monitor`.
271: Level: advanced
273: .seealso: [](ch:nep), `NEPMonitorSet()`, `NEPMonitorRegisterAll()`, `NEPMonitorSetFromOptions()`
274: @*/
275: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
276: {
277: char key[PETSC_MAX_PATH_LEN];
279: PetscFunctionBegin;
280: PetscCall(NEPInitializePackage());
281: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
282: PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
283: if (create) PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
284: if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*
289: NEPReset_Problem - Destroys the problem matrices.
290: */
291: PetscErrorCode NEPReset_Problem(NEP nep)
292: {
293: PetscInt i;
295: PetscFunctionBegin;
297: PetscCall(MatDestroy(&nep->function));
298: PetscCall(MatDestroy(&nep->function_pre));
299: PetscCall(MatDestroy(&nep->jacobian));
300: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
301: PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
302: for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
303: PetscCall(PetscFree(nep->f));
304: PetscCall(PetscFree(nep->nrma));
305: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
306: nep->nt = 0;
307: }
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
310: /*@
311: NEPReset - Resets the NEP context to the initial state (prior to setup)
312: and destroys any allocated Vecs and Mats.
314: Collective
316: Input Parameter:
317: . nep - the nonlinear eigensolver context
319: Level: advanced
321: .seealso: [](ch:nep), `NEPDestroy()`
322: @*/
323: PetscErrorCode NEPReset(NEP nep)
324: {
325: PetscFunctionBegin;
327: if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
328: PetscTryTypeMethod(nep,reset);
329: if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
330: PetscCall(NEPReset_Problem(nep));
331: PetscCall(BVDestroy(&nep->V));
332: PetscCall(BVDestroy(&nep->W));
333: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
334: PetscCall(MatDestroy(&nep->resolvent));
335: nep->nwork = 0;
336: nep->state = NEP_STATE_INITIAL;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: NEPDestroy - Destroys the NEP context.
343: Collective
345: Input Parameter:
346: . nep - the nonlinear eigensolver context
348: Level: beginner
350: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPSolve()`
351: @*/
352: PetscErrorCode NEPDestroy(NEP *nep)
353: {
354: PetscFunctionBegin;
355: if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
357: if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
358: PetscCall(NEPReset(*nep));
359: PetscTryTypeMethod(*nep,destroy);
360: if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
361: PetscCall(RGDestroy(&(*nep)->rg));
362: PetscCall(DSDestroy(&(*nep)->ds));
363: PetscCall(KSPDestroy(&(*nep)->refineksp));
364: PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
365: PetscCall(PetscFree((*nep)->sc));
366: /* just in case the initial vectors have not been used */
367: PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
368: if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
369: if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
370: PetscCall(NEPMonitorCancel(*nep));
371: PetscCall(PetscHeaderDestroy(nep));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /*@
376: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
378: Collective
380: Input Parameters:
381: + nep - the nonlinear eigensolver context
382: - bv - the basis vectors object
384: Note:
385: Use NEPGetBV() to retrieve the basis vectors context (for example,
386: to free it at the end of the computations).
388: Level: advanced
390: .seealso: [](ch:nep), `NEPGetBV()`
391: @*/
392: PetscErrorCode NEPSetBV(NEP nep,BV bv)
393: {
394: PetscFunctionBegin;
397: PetscCheckSameComm(nep,1,bv,2);
398: PetscCall(PetscObjectReference((PetscObject)bv));
399: PetscCall(BVDestroy(&nep->V));
400: nep->V = bv;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
406: eigensolver object.
408: Not Collective
410: Input Parameter:
411: . nep - the nonlinear eigensolver context
413: Output Parameter:
414: . bv - basis vectors context
416: Level: advanced
418: .seealso: [](ch:nep), `NEPSetBV()`
419: @*/
420: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
421: {
422: PetscFunctionBegin;
424: PetscAssertPointer(bv,2);
425: if (!nep->V) {
426: PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
427: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
428: PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
429: }
430: *bv = nep->V;
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: NEPSetRG - Associates a region object to the nonlinear eigensolver.
437: Collective
439: Input Parameters:
440: + nep - the nonlinear eigensolver context
441: - rg - the region object
443: Note:
444: Use NEPGetRG() to retrieve the region context (for example,
445: to free it at the end of the computations).
447: Level: advanced
449: .seealso: [](ch:nep), `NEPGetRG()`
450: @*/
451: PetscErrorCode NEPSetRG(NEP nep,RG rg)
452: {
453: PetscFunctionBegin;
455: if (rg) {
457: PetscCheckSameComm(nep,1,rg,2);
458: }
459: PetscCall(PetscObjectReference((PetscObject)rg));
460: PetscCall(RGDestroy(&nep->rg));
461: nep->rg = rg;
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: NEPGetRG - Obtain the region object associated to the
467: nonlinear eigensolver object.
469: Not Collective
471: Input Parameter:
472: . nep - the nonlinear eigensolver context
474: Output Parameter:
475: . rg - region context
477: Level: advanced
479: .seealso: [](ch:nep), `NEPSetRG()`
480: @*/
481: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
482: {
483: PetscFunctionBegin;
485: PetscAssertPointer(rg,2);
486: if (!nep->rg) {
487: PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
488: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
489: PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
490: }
491: *rg = nep->rg;
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
498: Collective
500: Input Parameters:
501: + nep - the nonlinear eigensolver context
502: - ds - the direct solver object
504: Note:
505: Use NEPGetDS() to retrieve the direct solver context (for example,
506: to free it at the end of the computations).
508: Level: advanced
510: .seealso: [](ch:nep), `NEPGetDS()`
511: @*/
512: PetscErrorCode NEPSetDS(NEP nep,DS ds)
513: {
514: PetscFunctionBegin;
517: PetscCheckSameComm(nep,1,ds,2);
518: PetscCall(PetscObjectReference((PetscObject)ds));
519: PetscCall(DSDestroy(&nep->ds));
520: nep->ds = ds;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: NEPGetDS - Obtain the direct solver object associated to the
526: nonlinear eigensolver object.
528: Not Collective
530: Input Parameter:
531: . nep - the nonlinear eigensolver context
533: Output Parameter:
534: . ds - direct solver context
536: Level: advanced
538: .seealso: [](ch:nep), `NEPSetDS()`
539: @*/
540: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
541: {
542: PetscFunctionBegin;
544: PetscAssertPointer(ds,2);
545: if (!nep->ds) {
546: PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
547: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
548: PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
549: }
550: *ds = nep->ds;
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: /*@
555: NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
556: object in the refinement phase.
558: Collective
560: Input Parameter:
561: . nep - the nonlinear eigensolver context
563: Output Parameter:
564: . ksp - ksp context
566: Level: advanced
568: .seealso: [](ch:nep), `NEPSetRefine()`
569: @*/
570: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
571: {
572: MPI_Comm comm;
574: PetscFunctionBegin;
576: PetscAssertPointer(ksp,2);
577: if (!nep->refineksp) {
578: if (nep->npart>1) {
579: /* Split in subcomunicators */
580: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
581: PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
582: PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
583: PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
584: } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
585: PetscCall(KSPCreate(comm,&nep->refineksp));
586: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
587: PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
588: PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
589: PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
590: PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
591: }
592: *ksp = nep->refineksp;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@
597: NEPSetTarget - Sets the value of the target.
599: Logically Collective
601: Input Parameters:
602: + nep - the nonlinear eigensolver context
603: - target - the value of the target
605: Options Database Key:
606: . -nep_target <scalar> - the value of the target
608: Notes:
609: The target is a scalar value used to determine the portion of the spectrum
610: of interest. It is used in combination with NEPSetWhichEigenpairs().
612: In the case of complex scalars, a complex value can be provided in the
613: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
614: -nep_target 1.0+2.0i
616: Level: intermediate
618: .seealso: [](ch:nep), `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
619: @*/
620: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
621: {
622: PetscFunctionBegin;
625: nep->target = target;
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /*@
630: NEPGetTarget - Gets the value of the target.
632: Not Collective
634: Input Parameter:
635: . nep - the nonlinear eigensolver context
637: Output Parameter:
638: . target - the value of the target
640: Note:
641: If the target was not set by the user, then zero is returned.
643: Level: intermediate
645: .seealso: [](ch:nep), `NEPSetTarget()`
646: @*/
647: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
648: {
649: PetscFunctionBegin;
651: PetscAssertPointer(target,2);
652: *target = nep->target;
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*@C
657: NEPSetFunction - Sets the function to compute the nonlinear Function T(lambda)
658: as well as the location to store the matrix.
660: Collective
662: Input Parameters:
663: + nep - the nonlinear eigensolver context
664: . F - Function matrix
665: . P - preconditioner matrix (usually same as F)
666: . fun - Function evaluation routine (if NULL then NEP retains any
667: previously set value), see NEPFunctionFn for the calling sequence
668: - ctx - [optional] user-defined context for private data for the Function
669: evaluation routine (may be NULL) (if NULL then NEP retains any
670: previously set value)
672: Level: beginner
674: .seealso: [](ch:nep), `NEPGetFunction()`, `NEPSetJacobian()`
675: @*/
676: PetscErrorCode NEPSetFunction(NEP nep,Mat F,Mat P,NEPFunctionFn *fun,void *ctx)
677: {
678: PetscFunctionBegin;
682: if (F) PetscCheckSameComm(nep,1,F,2);
683: if (P) PetscCheckSameComm(nep,1,P,3);
685: if (nep->state) PetscCall(NEPReset(nep));
686: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
688: if (fun) nep->computefunction = fun;
689: if (ctx) nep->functionctx = ctx;
690: if (F) {
691: PetscCall(PetscObjectReference((PetscObject)F));
692: PetscCall(MatDestroy(&nep->function));
693: nep->function = F;
694: }
695: if (P) {
696: PetscCall(PetscObjectReference((PetscObject)P));
697: PetscCall(MatDestroy(&nep->function_pre));
698: nep->function_pre = P;
699: }
700: nep->fui = NEP_USER_INTERFACE_CALLBACK;
701: nep->state = NEP_STATE_INITIAL;
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: /*@C
706: NEPGetFunction - Returns the Function matrix and optionally the user
707: provided context for evaluating the Function.
709: Collective
711: Input Parameter:
712: . nep - the nonlinear eigensolver context
714: Output Parameters:
715: + F - location to stash Function matrix (or NULL)
716: . P - location to stash preconditioner matrix (or NULL)
717: . fun - location to put Function function (or NULL)
718: - ctx - location to stash Function context (or NULL)
720: Level: advanced
722: .seealso: [](ch:nep), `NEPSetFunction()`
723: @*/
724: PetscErrorCode NEPGetFunction(NEP nep,Mat *F,Mat *P,NEPFunctionFn **fun,void **ctx)
725: {
726: PetscFunctionBegin;
728: NEPCheckCallback(nep,1);
729: if (F) *F = nep->function;
730: if (P) *P = nep->function_pre;
731: if (fun) *fun = nep->computefunction;
732: if (ctx) *ctx = nep->functionctx;
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*@C
737: NEPSetJacobian - Sets the function to compute the Jacobian T'(lambda) as well
738: as the location to store the matrix.
740: Collective
742: Input Parameters:
743: + nep - the nonlinear eigensolver context
744: . J - Jacobian matrix
745: . jac - Jacobian evaluation routine (if NULL then NEP retains any
746: previously set value), see NEPJacobianFn for the calling sequence
747: - ctx - [optional] user-defined context for private data for the Jacobian
748: evaluation routine (may be NULL) (if NULL then NEP retains any
749: previously set value)
751: Level: beginner
753: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetJacobian()`
754: @*/
755: PetscErrorCode NEPSetJacobian(NEP nep,Mat J,NEPJacobianFn *jac,void *ctx)
756: {
757: PetscFunctionBegin;
760: if (J) PetscCheckSameComm(nep,1,J,2);
762: if (nep->state) PetscCall(NEPReset(nep));
763: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
765: if (jac) nep->computejacobian = jac;
766: if (ctx) nep->jacobianctx = ctx;
767: if (J) {
768: PetscCall(PetscObjectReference((PetscObject)J));
769: PetscCall(MatDestroy(&nep->jacobian));
770: nep->jacobian = J;
771: }
772: nep->fui = NEP_USER_INTERFACE_CALLBACK;
773: nep->state = NEP_STATE_INITIAL;
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: /*@C
778: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
779: provided routine and context for evaluating the Jacobian.
781: Collective
783: Input Parameter:
784: . nep - the nonlinear eigensolver context
786: Output Parameters:
787: + J - location to stash Jacobian matrix (or NULL)
788: . jac - location to put Jacobian function (or NULL)
789: - ctx - location to stash Jacobian context (or NULL)
791: Level: advanced
793: .seealso: [](ch:nep), `NEPSetJacobian()`
794: @*/
795: PetscErrorCode NEPGetJacobian(NEP nep,Mat *J,NEPJacobianFn **jac,void **ctx)
796: {
797: PetscFunctionBegin;
799: NEPCheckCallback(nep,1);
800: if (J) *J = nep->jacobian;
801: if (jac) *jac = nep->computejacobian;
802: if (ctx) *ctx = nep->jacobianctx;
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: /*@
807: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
808: in split form.
810: Collective
812: Input Parameters:
813: + nep - the nonlinear eigensolver context
814: . nt - number of terms in the split form
815: . A - array of matrices
816: . f - array of functions
817: - str - structure flag for matrices
819: Notes:
820: The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
821: for i=1,...,n. The derivative T'(lambda) can be obtained using the
822: derivatives of f_i.
824: The structure flag provides information about A_i's nonzero pattern
825: (see MatStructure enum). If all matrices have the same pattern, then
826: use SAME_NONZERO_PATTERN. If the patterns are different but contained
827: in the pattern of the first one, then use SUBSET_NONZERO_PATTERN. If
828: patterns are known to be different, use DIFFERENT_NONZERO_PATTERN.
829: If set to UNKNOWN_NONZERO_PATTERN, the patterns will be compared to
830: determine if they are equal.
832: This function must be called before NEPSetUp(). If it is called again
833: after NEPSetUp() then the NEP object is reset.
835: Level: beginner
837: .seealso: [](ch:nep), `NEPGetSplitOperatorTerm()`, `NEPGetSplitOperatorInfo()`, `NEPSetSplitPreconditioner()`
838: @*/
839: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
840: {
841: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
843: PetscFunctionBegin;
846: PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
847: PetscAssertPointer(A,3);
848: PetscAssertPointer(f,4);
851: for (i=0;i<nt;i++) {
853: PetscCheckSameComm(nep,1,A[i],3);
855: PetscCheckSameComm(nep,1,f[i],4);
856: PetscCall(MatGetSize(A[i],&m,&n));
857: PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
858: PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
859: PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
860: if (!i) { m0 = m; mloc0 = mloc; }
861: PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
862: PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
863: PetscCall(PetscObjectReference((PetscObject)A[i]));
864: PetscCall(PetscObjectReference((PetscObject)f[i]));
865: }
867: if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
868: else PetscCall(NEPReset_Problem(nep));
870: /* allocate space and copy matrices and functions */
871: PetscCall(PetscMalloc1(nt,&nep->A));
872: for (i=0;i<nt;i++) nep->A[i] = A[i];
873: PetscCall(PetscMalloc1(nt,&nep->f));
874: for (i=0;i<nt;i++) nep->f[i] = f[i];
875: PetscCall(PetscCalloc1(nt,&nep->nrma));
876: nep->nt = nt;
877: nep->mstr = str;
878: nep->fui = NEP_USER_INTERFACE_SPLIT;
879: nep->state = NEP_STATE_INITIAL;
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: /*@
884: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
885: the nonlinear operator in split form.
887: Collective
889: Input Parameters:
890: + nep - the nonlinear eigensolver context
891: - k - the index of the requested term (starting in 0)
893: Output Parameters:
894: + A - the matrix of the requested term
895: - f - the function of the requested term
897: Level: intermediate
899: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorInfo()`
900: @*/
901: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
902: {
903: PetscFunctionBegin;
906: NEPCheckSplit(nep,1);
907: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
908: if (A) *A = nep->A[k];
909: if (f) *f = nep->f[k];
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@
914: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
915: the nonlinear operator, as well as the structure flag for matrices.
917: Not Collective
919: Input Parameter:
920: . nep - the nonlinear eigensolver context
922: Output Parameters:
923: + n - the number of terms passed in NEPSetSplitOperator()
924: - str - the matrix structure flag passed in NEPSetSplitOperator()
926: Level: intermediate
928: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorTerm()`
929: @*/
930: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
931: {
932: PetscFunctionBegin;
934: NEPCheckSplit(nep,1);
935: if (n) *n = nep->nt;
936: if (str) *str = nep->mstr;
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: /*@
941: NEPSetSplitPreconditioner - Sets an operator in split form from which
942: to build the preconditioner to be used when solving the nonlinear
943: eigenvalue problem in split form.
945: Collective
947: Input Parameters:
948: + nep - the nonlinear eigensolver context
949: . ntp - number of terms in the split preconditioner
950: . P - array of matrices
951: - strp - structure flag for matrices
953: Notes:
954: The matrix for the preconditioner is expressed as P(lambda) =
955: sum_i P_i*f_i(lambda), for i=1,...,n, where the f_i functions
956: are the same as in NEPSetSplitOperator(). It is not necessary to call
957: this function. If it is not invoked, then the preconditioner is
958: built from T(lambda), i.e., both matrices and functions passed in
959: NEPSetSplitOperator().
961: The structure flag provides information about P_i's nonzero pattern
962: in the same way as in NEPSetSplitOperator().
964: If the functions defining the preconditioner operator were different
965: from the ones given in NEPSetSplitOperator(), then the split form
966: cannot be used. Use the callback interface instead.
968: Use ntp=0 to reset a previously set split preconditioner.
970: Level: advanced
972: .seealso: [](ch:nep), `NEPGetSplitPreconditionerTerm()`, `NEPGetSplitPreconditionerInfo()`, `NEPSetSplitOperator()`
973: @*/
974: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
975: {
976: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
978: PetscFunctionBegin;
981: PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
982: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
983: PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
984: if (ntp) PetscAssertPointer(P,3);
987: for (i=0;i<ntp;i++) {
989: PetscCheckSameComm(nep,1,P[i],3);
990: PetscCall(MatGetSize(P[i],&m,&n));
991: PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
992: PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
993: PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
994: if (!i) { m0 = m; mloc0 = mloc; }
995: PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
996: PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
997: PetscCall(PetscObjectReference((PetscObject)P[i]));
998: }
1000: PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
1001: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
1003: /* allocate space and copy matrices */
1004: if (ntp) {
1005: PetscCall(PetscMalloc1(ntp,&nep->P));
1006: for (i=0;i<ntp;i++) nep->P[i] = P[i];
1007: }
1008: nep->mstrp = strp;
1009: nep->state = NEP_STATE_INITIAL;
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*@
1014: NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1015: the split preconditioner.
1017: Not Collective
1019: Input Parameters:
1020: + nep - the nonlinear eigensolver context
1021: - k - the index of the requested term (starting in 0)
1023: Output Parameter:
1024: . P - the matrix of the requested term
1026: Level: advanced
1028: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerInfo()`
1029: @*/
1030: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1031: {
1032: PetscFunctionBegin;
1035: PetscAssertPointer(P,3);
1036: NEPCheckSplit(nep,1);
1037: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1038: PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1039: *P = nep->P[k];
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: /*@
1044: NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1045: preconditioner, as well as the structure flag for matrices.
1047: Not Collective
1049: Input Parameter:
1050: . nep - the nonlinear eigensolver context
1052: Output Parameters:
1053: + n - the number of terms passed in NEPSetSplitPreconditioner()
1054: - strp - the matrix structure flag passed in NEPSetSplitPreconditioner()
1056: Level: advanced
1058: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerTerm()`
1059: @*/
1060: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1061: {
1062: PetscFunctionBegin;
1064: NEPCheckSplit(nep,1);
1065: if (n) *n = nep->P? nep->nt: 0;
1066: if (strp) *strp = nep->mstrp;
1067: PetscFunctionReturn(PETSC_SUCCESS);
1068: }