Actual source code: nepbasic.c
slepc-main 2024-11-22
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 - Adds NEP monitor routine.
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
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: Example Usage:
252: .vb
253: NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
254: .ve
256: Then, your monitor can be chosen with the procedural interface via
257: $ NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL)
258: or at runtime via the option
259: $ -nep_monitor_my_monitor
261: Level: advanced
263: .seealso: NEPMonitorRegisterAll()
264: @*/
265: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(NEP,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
266: {
267: char key[PETSC_MAX_PATH_LEN];
269: PetscFunctionBegin;
270: PetscCall(NEPInitializePackage());
271: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
272: PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
273: if (create) PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
274: if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /*
279: NEPReset_Problem - Destroys the problem matrices.
280: */
281: PetscErrorCode NEPReset_Problem(NEP nep)
282: {
283: PetscInt i;
285: PetscFunctionBegin;
287: PetscCall(MatDestroy(&nep->function));
288: PetscCall(MatDestroy(&nep->function_pre));
289: PetscCall(MatDestroy(&nep->jacobian));
290: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
291: PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
292: for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
293: PetscCall(PetscFree(nep->f));
294: PetscCall(PetscFree(nep->nrma));
295: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
296: nep->nt = 0;
297: }
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
300: /*@
301: NEPReset - Resets the NEP context to the initial state (prior to setup)
302: and destroys any allocated Vecs and Mats.
304: Collective
306: Input Parameter:
307: . nep - eigensolver context obtained from NEPCreate()
309: Level: advanced
311: .seealso: NEPDestroy()
312: @*/
313: PetscErrorCode NEPReset(NEP nep)
314: {
315: PetscFunctionBegin;
317: if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
318: PetscTryTypeMethod(nep,reset);
319: if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
320: PetscCall(NEPReset_Problem(nep));
321: PetscCall(BVDestroy(&nep->V));
322: PetscCall(BVDestroy(&nep->W));
323: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
324: PetscCall(MatDestroy(&nep->resolvent));
325: nep->nwork = 0;
326: nep->state = NEP_STATE_INITIAL;
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*@
331: NEPDestroy - Destroys the NEP context.
333: Collective
335: Input Parameter:
336: . nep - eigensolver context obtained from NEPCreate()
338: Level: beginner
340: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
341: @*/
342: PetscErrorCode NEPDestroy(NEP *nep)
343: {
344: PetscFunctionBegin;
345: if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
347: if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
348: PetscCall(NEPReset(*nep));
349: PetscTryTypeMethod(*nep,destroy);
350: if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
351: PetscCall(RGDestroy(&(*nep)->rg));
352: PetscCall(DSDestroy(&(*nep)->ds));
353: PetscCall(KSPDestroy(&(*nep)->refineksp));
354: PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
355: PetscCall(PetscFree((*nep)->sc));
356: /* just in case the initial vectors have not been used */
357: PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
358: if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
359: if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
360: PetscCall(NEPMonitorCancel(*nep));
361: PetscCall(PetscHeaderDestroy(nep));
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
368: Collective
370: Input Parameters:
371: + nep - eigensolver context obtained from NEPCreate()
372: - bv - the basis vectors object
374: Note:
375: Use NEPGetBV() to retrieve the basis vectors context (for example,
376: to free it at the end of the computations).
378: Level: advanced
380: .seealso: NEPGetBV()
381: @*/
382: PetscErrorCode NEPSetBV(NEP nep,BV bv)
383: {
384: PetscFunctionBegin;
387: PetscCheckSameComm(nep,1,bv,2);
388: PetscCall(PetscObjectReference((PetscObject)bv));
389: PetscCall(BVDestroy(&nep->V));
390: nep->V = bv;
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@
395: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
396: eigensolver object.
398: Not Collective
400: Input Parameters:
401: . nep - eigensolver context obtained from NEPCreate()
403: Output Parameter:
404: . bv - basis vectors context
406: Level: advanced
408: .seealso: NEPSetBV()
409: @*/
410: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
411: {
412: PetscFunctionBegin;
414: PetscAssertPointer(bv,2);
415: if (!nep->V) {
416: PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
417: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
418: PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
419: }
420: *bv = nep->V;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: NEPSetRG - Associates a region object to the nonlinear eigensolver.
427: Collective
429: Input Parameters:
430: + nep - eigensolver context obtained from NEPCreate()
431: - rg - the region object
433: Note:
434: Use NEPGetRG() to retrieve the region context (for example,
435: to free it at the end of the computations).
437: Level: advanced
439: .seealso: NEPGetRG()
440: @*/
441: PetscErrorCode NEPSetRG(NEP nep,RG rg)
442: {
443: PetscFunctionBegin;
445: if (rg) {
447: PetscCheckSameComm(nep,1,rg,2);
448: }
449: PetscCall(PetscObjectReference((PetscObject)rg));
450: PetscCall(RGDestroy(&nep->rg));
451: nep->rg = rg;
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: /*@
456: NEPGetRG - Obtain the region object associated to the
457: nonlinear eigensolver object.
459: Not Collective
461: Input Parameters:
462: . nep - eigensolver context obtained from NEPCreate()
464: Output Parameter:
465: . rg - region context
467: Level: advanced
469: .seealso: NEPSetRG()
470: @*/
471: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
472: {
473: PetscFunctionBegin;
475: PetscAssertPointer(rg,2);
476: if (!nep->rg) {
477: PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
478: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
479: PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
480: }
481: *rg = nep->rg;
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: /*@
486: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
488: Collective
490: Input Parameters:
491: + nep - eigensolver context obtained from NEPCreate()
492: - ds - the direct solver object
494: Note:
495: Use NEPGetDS() to retrieve the direct solver context (for example,
496: to free it at the end of the computations).
498: Level: advanced
500: .seealso: NEPGetDS()
501: @*/
502: PetscErrorCode NEPSetDS(NEP nep,DS ds)
503: {
504: PetscFunctionBegin;
507: PetscCheckSameComm(nep,1,ds,2);
508: PetscCall(PetscObjectReference((PetscObject)ds));
509: PetscCall(DSDestroy(&nep->ds));
510: nep->ds = ds;
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*@
515: NEPGetDS - Obtain the direct solver object associated to the
516: nonlinear eigensolver object.
518: Not Collective
520: Input Parameters:
521: . nep - eigensolver context obtained from NEPCreate()
523: Output Parameter:
524: . ds - direct solver context
526: Level: advanced
528: .seealso: NEPSetDS()
529: @*/
530: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
531: {
532: PetscFunctionBegin;
534: PetscAssertPointer(ds,2);
535: if (!nep->ds) {
536: PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
537: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
538: PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
539: }
540: *ds = nep->ds;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
546: object in the refinement phase.
548: Collective
550: Input Parameters:
551: . nep - eigensolver context obtained from NEPCreate()
553: Output Parameter:
554: . ksp - ksp context
556: Level: advanced
558: .seealso: NEPSetRefine()
559: @*/
560: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
561: {
562: MPI_Comm comm;
564: PetscFunctionBegin;
566: PetscAssertPointer(ksp,2);
567: if (!nep->refineksp) {
568: if (nep->npart>1) {
569: /* Split in subcomunicators */
570: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
571: PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
572: PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
573: PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
574: } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
575: PetscCall(KSPCreate(comm,&nep->refineksp));
576: PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
577: PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
578: PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
579: PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
580: PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
581: }
582: *ksp = nep->refineksp;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: NEPSetTarget - Sets the value of the target.
589: Logically Collective
591: Input Parameters:
592: + nep - eigensolver context
593: - target - the value of the target
595: Options Database Key:
596: . -nep_target <scalar> - the value of the target
598: Notes:
599: The target is a scalar value used to determine the portion of the spectrum
600: of interest. It is used in combination with NEPSetWhichEigenpairs().
602: In the case of complex scalars, a complex value can be provided in the
603: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
604: -nep_target 1.0+2.0i
606: Level: intermediate
608: .seealso: NEPGetTarget(), NEPSetWhichEigenpairs()
609: @*/
610: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
611: {
612: PetscFunctionBegin;
615: nep->target = target;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@
620: NEPGetTarget - Gets the value of the target.
622: Not Collective
624: Input Parameter:
625: . nep - eigensolver context
627: Output Parameter:
628: . target - the value of the target
630: Note:
631: If the target was not set by the user, then zero is returned.
633: Level: intermediate
635: .seealso: NEPSetTarget()
636: @*/
637: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
638: {
639: PetscFunctionBegin;
641: PetscAssertPointer(target,2);
642: *target = nep->target;
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: /*@C
647: NEPSetFunction - Sets the function to compute the nonlinear Function T(lambda)
648: as well as the location to store the matrix.
650: Collective
652: Input Parameters:
653: + nep - the NEP context
654: . A - Function matrix
655: . B - preconditioner matrix (usually same as A)
656: . fun - Function evaluation routine (if NULL then NEP retains any
657: previously set value), see NEPFunctionFn for the calling sequence
658: - ctx - [optional] user-defined context for private data for the Function
659: evaluation routine (may be NULL) (if NULL then NEP retains any
660: previously set value)
662: Level: beginner
664: .seealso: NEPGetFunction(), NEPSetJacobian()
665: @*/
666: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,NEPFunctionFn *fun,void *ctx)
667: {
668: PetscFunctionBegin;
672: if (A) PetscCheckSameComm(nep,1,A,2);
673: if (B) PetscCheckSameComm(nep,1,B,3);
675: if (nep->state) PetscCall(NEPReset(nep));
676: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
678: if (fun) nep->computefunction = fun;
679: if (ctx) nep->functionctx = ctx;
680: if (A) {
681: PetscCall(PetscObjectReference((PetscObject)A));
682: PetscCall(MatDestroy(&nep->function));
683: nep->function = A;
684: }
685: if (B) {
686: PetscCall(PetscObjectReference((PetscObject)B));
687: PetscCall(MatDestroy(&nep->function_pre));
688: nep->function_pre = B;
689: }
690: nep->fui = NEP_USER_INTERFACE_CALLBACK;
691: nep->state = NEP_STATE_INITIAL;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@C
696: NEPGetFunction - Returns the Function matrix and optionally the user
697: provided context for evaluating the Function.
699: Not Collective
701: Input Parameter:
702: . nep - the nonlinear eigensolver context
704: Output Parameters:
705: + A - location to stash Function matrix (or NULL)
706: . B - location to stash preconditioner matrix (or NULL)
707: . fun - location to put Function function (or NULL)
708: - ctx - location to stash Function context (or NULL)
710: Level: advanced
712: .seealso: NEPSetFunction()
713: @*/
714: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,NEPFunctionFn **fun,void **ctx)
715: {
716: PetscFunctionBegin;
718: NEPCheckCallback(nep,1);
719: if (A) *A = nep->function;
720: if (B) *B = nep->function_pre;
721: if (fun) *fun = nep->computefunction;
722: if (ctx) *ctx = nep->functionctx;
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@C
727: NEPSetJacobian - Sets the function to compute the Jacobian T'(lambda) as well
728: as the location to store the matrix.
730: Collective
732: Input Parameters:
733: + nep - the NEP context
734: . A - Jacobian matrix
735: . jac - Jacobian evaluation routine (if NULL then NEP retains any
736: previously set value), see NEPJacobianFn for the calling sequence
737: - ctx - [optional] user-defined context for private data for the Jacobian
738: evaluation routine (may be NULL) (if NULL then NEP retains any
739: previously set value)
741: Level: beginner
743: .seealso: NEPSetFunction(), NEPGetJacobian()
744: @*/
745: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,NEPJacobianFn *jac,void *ctx)
746: {
747: PetscFunctionBegin;
750: if (A) PetscCheckSameComm(nep,1,A,2);
752: if (nep->state) PetscCall(NEPReset(nep));
753: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
755: if (jac) nep->computejacobian = jac;
756: if (ctx) nep->jacobianctx = ctx;
757: if (A) {
758: PetscCall(PetscObjectReference((PetscObject)A));
759: PetscCall(MatDestroy(&nep->jacobian));
760: nep->jacobian = A;
761: }
762: nep->fui = NEP_USER_INTERFACE_CALLBACK;
763: nep->state = NEP_STATE_INITIAL;
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: /*@C
768: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
769: provided routine and context for evaluating the Jacobian.
771: Not Collective
773: Input Parameter:
774: . nep - the nonlinear eigensolver context
776: Output Parameters:
777: + A - location to stash Jacobian matrix (or NULL)
778: . jac - location to put Jacobian function (or NULL)
779: - ctx - location to stash Jacobian context (or NULL)
781: Level: advanced
783: .seealso: NEPSetJacobian()
784: @*/
785: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,NEPJacobianFn **jac,void **ctx)
786: {
787: PetscFunctionBegin;
789: NEPCheckCallback(nep,1);
790: if (A) *A = nep->jacobian;
791: if (jac) *jac = nep->computejacobian;
792: if (ctx) *ctx = nep->jacobianctx;
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }
796: /*@
797: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
798: in split form.
800: Collective
802: Input Parameters:
803: + nep - the nonlinear eigensolver context
804: . nt - number of terms in the split form
805: . A - array of matrices
806: . f - array of functions
807: - str - structure flag for matrices
809: Notes:
810: The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
811: for i=1,...,n. The derivative T'(lambda) can be obtained using the
812: derivatives of f_i.
814: The structure flag provides information about A_i's nonzero pattern
815: (see MatStructure enum). If all matrices have the same pattern, then
816: use SAME_NONZERO_PATTERN. If the patterns are different but contained
817: in the pattern of the first one, then use SUBSET_NONZERO_PATTERN. If
818: patterns are known to be different, use DIFFERENT_NONZERO_PATTERN.
819: If set to UNKNOWN_NONZERO_PATTERN, the patterns will be compared to
820: determine if they are equal.
822: This function must be called before NEPSetUp(). If it is called again
823: after NEPSetUp() then the NEP object is reset.
825: Level: beginner
827: .seealso: NEPGetSplitOperatorTerm(), NEPGetSplitOperatorInfo(), NEPSetSplitPreconditioner()
828: @*/
829: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
830: {
831: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
833: PetscFunctionBegin;
836: PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
837: PetscAssertPointer(A,3);
838: PetscAssertPointer(f,4);
841: for (i=0;i<nt;i++) {
843: PetscCheckSameComm(nep,1,A[i],3);
845: PetscCheckSameComm(nep,1,f[i],4);
846: PetscCall(MatGetSize(A[i],&m,&n));
847: PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
848: 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);
849: 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);
850: if (!i) { m0 = m; mloc0 = mloc; }
851: 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);
852: 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);
853: PetscCall(PetscObjectReference((PetscObject)A[i]));
854: PetscCall(PetscObjectReference((PetscObject)f[i]));
855: }
857: if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
858: else PetscCall(NEPReset_Problem(nep));
860: /* allocate space and copy matrices and functions */
861: PetscCall(PetscMalloc1(nt,&nep->A));
862: for (i=0;i<nt;i++) nep->A[i] = A[i];
863: PetscCall(PetscMalloc1(nt,&nep->f));
864: for (i=0;i<nt;i++) nep->f[i] = f[i];
865: PetscCall(PetscCalloc1(nt,&nep->nrma));
866: nep->nt = nt;
867: nep->mstr = str;
868: nep->fui = NEP_USER_INTERFACE_SPLIT;
869: nep->state = NEP_STATE_INITIAL;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
875: the nonlinear operator in split form.
877: Not Collective
879: Input Parameters:
880: + nep - the nonlinear eigensolver context
881: - k - the index of the requested term (starting in 0)
883: Output Parameters:
884: + A - the matrix of the requested term
885: - f - the function of the requested term
887: Level: intermediate
889: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
890: @*/
891: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
892: {
893: PetscFunctionBegin;
896: NEPCheckSplit(nep,1);
897: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
898: if (A) *A = nep->A[k];
899: if (f) *f = nep->f[k];
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: /*@
904: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
905: the nonlinear operator, as well as the structure flag for matrices.
907: Not Collective
909: Input Parameter:
910: . nep - the nonlinear eigensolver context
912: Output Parameters:
913: + n - the number of terms passed in NEPSetSplitOperator()
914: - str - the matrix structure flag passed in NEPSetSplitOperator()
916: Level: intermediate
918: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
919: @*/
920: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
921: {
922: PetscFunctionBegin;
924: NEPCheckSplit(nep,1);
925: if (n) *n = nep->nt;
926: if (str) *str = nep->mstr;
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: /*@
931: NEPSetSplitPreconditioner - Sets an operator in split form from which
932: to build the preconditioner to be used when solving the nonlinear
933: eigenvalue problem in split form.
935: Collective
937: Input Parameters:
938: + nep - the nonlinear eigensolver context
939: . ntp - number of terms in the split preconditioner
940: . P - array of matrices
941: - strp - structure flag for matrices
943: Notes:
944: The matrix for the preconditioner is expressed as P(lambda) =
945: sum_i P_i*f_i(lambda), for i=1,...,n, where the f_i functions
946: are the same as in NEPSetSplitOperator(). It is not necessary to call
947: this function. If it is not invoked, then the preconditioner is
948: built from T(lambda), i.e., both matrices and functions passed in
949: NEPSetSplitOperator().
951: The structure flag provides information about P_i's nonzero pattern
952: in the same way as in NEPSetSplitOperator().
954: If the functions defining the preconditioner operator were different
955: from the ones given in NEPSetSplitOperator(), then the split form
956: cannot be used. Use the callback interface instead.
958: Use ntp=0 to reset a previously set split preconditioner.
960: Level: advanced
962: .seealso: NEPGetSplitPreconditionerTerm(), NEPGetSplitPreconditionerInfo(), NEPSetSplitOperator()
963: @*/
964: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
965: {
966: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
968: PetscFunctionBegin;
971: PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
972: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
973: PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
974: if (ntp) PetscAssertPointer(P,3);
977: for (i=0;i<ntp;i++) {
979: PetscCheckSameComm(nep,1,P[i],3);
980: PetscCall(MatGetSize(P[i],&m,&n));
981: PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
982: 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);
983: 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);
984: if (!i) { m0 = m; mloc0 = mloc; }
985: 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);
986: 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);
987: PetscCall(PetscObjectReference((PetscObject)P[i]));
988: }
990: PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
991: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
993: /* allocate space and copy matrices */
994: if (ntp) {
995: PetscCall(PetscMalloc1(ntp,&nep->P));
996: for (i=0;i<ntp;i++) nep->P[i] = P[i];
997: }
998: nep->mstrp = strp;
999: nep->state = NEP_STATE_INITIAL;
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1005: the split preconditioner.
1007: Not Collective
1009: Input Parameters:
1010: + nep - the nonlinear eigensolver context
1011: - k - the index of the requested term (starting in 0)
1013: Output Parameter:
1014: . P - the matrix of the requested term
1016: Level: advanced
1018: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerInfo()
1019: @*/
1020: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1021: {
1022: PetscFunctionBegin;
1025: PetscAssertPointer(P,3);
1026: NEPCheckSplit(nep,1);
1027: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1028: PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1029: *P = nep->P[k];
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: /*@
1034: NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1035: preconditioner, as well as the structure flag for matrices.
1037: Not Collective
1039: Input Parameter:
1040: . nep - the nonlinear eigensolver context
1042: Output Parameters:
1043: + n - the number of terms passed in NEPSetSplitPreconditioner()
1044: - strp - the matrix structure flag passed in NEPSetSplitPreconditioner()
1046: Level: advanced
1048: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerTerm()
1049: @*/
1050: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1051: {
1052: PetscFunctionBegin;
1054: NEPCheckSplit(nep,1);
1055: if (n) *n = nep->P? nep->nt: 0;
1056: if (strp) *strp = nep->mstrp;
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }