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