Actual source code: nepbasic.c
slepc-3.21.0 2024-03-30
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: *outnep = NULL;
52: PetscCall(NEPInitializePackage());
53: PetscCall(SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView));
55: nep->max_it = PETSC_DEFAULT;
56: nep->nev = 1;
57: nep->ncv = PETSC_DEFAULT;
58: nep->mpd = PETSC_DEFAULT;
59: nep->nini = 0;
60: nep->target = 0.0;
61: nep->tol = PETSC_DEFAULT;
62: nep->conv = NEP_CONV_REL;
63: nep->stop = NEP_STOP_BASIC;
64: nep->which = (NEPWhich)0;
65: nep->problem_type = (NEPProblemType)0;
66: nep->refine = NEP_REFINE_NONE;
67: nep->npart = 1;
68: nep->rtol = PETSC_DEFAULT;
69: nep->rits = PETSC_DEFAULT;
70: nep->scheme = (NEPRefineScheme)0;
71: nep->trackall = PETSC_FALSE;
72: nep->twosided = PETSC_FALSE;
74: nep->computefunction = NULL;
75: nep->computejacobian = NULL;
76: nep->functionctx = NULL;
77: nep->jacobianctx = NULL;
78: nep->converged = NEPConvergedRelative;
79: nep->convergeduser = NULL;
80: nep->convergeddestroy= NULL;
81: nep->stopping = NEPStoppingBasic;
82: nep->stoppinguser = NULL;
83: nep->stoppingdestroy = NULL;
84: nep->convergedctx = NULL;
85: nep->stoppingctx = NULL;
86: nep->numbermonitors = 0;
88: nep->ds = NULL;
89: nep->V = NULL;
90: nep->W = NULL;
91: nep->rg = NULL;
92: nep->function = NULL;
93: nep->function_pre = NULL;
94: nep->jacobian = NULL;
95: nep->A = NULL;
96: nep->f = NULL;
97: nep->nt = 0;
98: nep->mstr = UNKNOWN_NONZERO_PATTERN;
99: nep->P = NULL;
100: nep->mstrp = UNKNOWN_NONZERO_PATTERN;
101: nep->IS = NULL;
102: nep->eigr = NULL;
103: nep->eigi = NULL;
104: nep->errest = NULL;
105: nep->perm = NULL;
106: nep->nwork = 0;
107: nep->work = NULL;
108: nep->data = NULL;
110: nep->state = NEP_STATE_INITIAL;
111: nep->nconv = 0;
112: nep->its = 0;
113: nep->n = 0;
114: nep->nloc = 0;
115: nep->nrma = NULL;
116: nep->fui = (NEPUserInterface)0;
117: nep->useds = PETSC_FALSE;
118: nep->resolvent = NULL;
119: nep->reason = NEP_CONVERGED_ITERATING;
121: PetscCall(PetscNew(&nep->sc));
122: *outnep = nep;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@C
127: NEPSetType - Selects the particular solver to be used in the NEP object.
129: Logically Collective
131: Input Parameters:
132: + nep - the nonlinear eigensolver context
133: - type - a known method
135: Options Database Key:
136: . -nep_type <method> - Sets the method; use -help for a list
137: of available methods
139: Notes:
140: See "slepc/include/slepcnep.h" for available methods.
142: Normally, it is best to use the NEPSetFromOptions() command and
143: then set the NEP type from the options database rather than by using
144: this routine. Using the options database provides the user with
145: maximum flexibility in evaluating the different available methods.
146: The NEPSetType() routine is provided for those situations where it
147: is necessary to set the iterative solver independently of the command
148: line or options database.
150: Level: intermediate
152: .seealso: NEPType
153: @*/
154: PetscErrorCode NEPSetType(NEP nep,NEPType type)
155: {
156: PetscErrorCode (*r)(NEP);
157: PetscBool match;
159: PetscFunctionBegin;
161: PetscAssertPointer(type,2);
163: PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
164: if (match) PetscFunctionReturn(PETSC_SUCCESS);
166: PetscCall(PetscFunctionListFind(NEPList,type,&r));
167: PetscCheck(r,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);
169: PetscTryTypeMethod(nep,destroy);
170: PetscCall(PetscMemzero(nep->ops,sizeof(struct _NEPOps)));
172: nep->state = NEP_STATE_INITIAL;
173: PetscCall(PetscObjectChangeTypeName((PetscObject)nep,type));
174: PetscCall((*r)(nep));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@C
179: NEPGetType - Gets the NEP type as a string from the NEP object.
181: Not Collective
183: Input Parameter:
184: . nep - the eigensolver context
186: Output Parameter:
187: . type - name of NEP method
189: Level: intermediate
191: .seealso: NEPSetType()
192: @*/
193: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
194: {
195: PetscFunctionBegin;
197: PetscAssertPointer(type,2);
198: *type = ((PetscObject)nep)->type_name;
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: /*@C
203: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
205: Not Collective
207: Input Parameters:
208: + name - name of a new user-defined solver
209: - function - routine to create the solver context
211: Notes:
212: NEPRegister() may be called multiple times to add several user-defined solvers.
214: Example Usage:
215: .vb
216: NEPRegister("my_solver",MySolverCreate);
217: .ve
219: Then, your solver can be chosen with the procedural interface via
220: $ NEPSetType(nep,"my_solver")
221: or at runtime via the option
222: $ -nep_type my_solver
224: Level: advanced
226: .seealso: NEPRegisterAll()
227: @*/
228: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
229: {
230: PetscFunctionBegin;
231: PetscCall(NEPInitializePackage());
232: PetscCall(PetscFunctionListAdd(&NEPList,name,function));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: /*@C
237: NEPMonitorRegister - Adds NEP monitor routine.
239: Not Collective
241: Input Parameters:
242: + name - name of a new monitor routine
243: . vtype - a PetscViewerType for the output
244: . format - a PetscViewerFormat for the output
245: . monitor - monitor routine
246: . create - creation routine, or NULL
247: - destroy - destruction routine, or NULL
249: Notes:
250: NEPMonitorRegister() may be called multiple times to add several user-defined monitors.
252: Example Usage:
253: .vb
254: NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
255: .ve
257: Then, your monitor can be chosen with the procedural interface via
258: $ NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL)
259: or at runtime via the option
260: $ -nep_monitor_my_monitor
262: Level: advanced
264: .seealso: NEPMonitorRegisterAll()
265: @*/
266: 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**))
267: {
268: char key[PETSC_MAX_PATH_LEN];
270: PetscFunctionBegin;
271: PetscCall(NEPInitializePackage());
272: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
273: PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
274: if (create) PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
275: if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*
280: NEPReset_Problem - Destroys the problem matrices.
281: */
282: PetscErrorCode NEPReset_Problem(NEP nep)
283: {
284: PetscInt i;
286: PetscFunctionBegin;
288: PetscCall(MatDestroy(&nep->function));
289: PetscCall(MatDestroy(&nep->function_pre));
290: PetscCall(MatDestroy(&nep->jacobian));
291: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
292: PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
293: for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
294: PetscCall(PetscFree(nep->f));
295: PetscCall(PetscFree(nep->nrma));
296: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
297: nep->nt = 0;
298: }
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
301: /*@
302: NEPReset - Resets the NEP context to the initial state (prior to setup)
303: and destroys any allocated Vecs and Mats.
305: Collective
307: Input Parameter:
308: . nep - eigensolver context obtained from NEPCreate()
310: Level: advanced
312: .seealso: NEPDestroy()
313: @*/
314: PetscErrorCode NEPReset(NEP nep)
315: {
316: PetscFunctionBegin;
318: if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
319: PetscTryTypeMethod(nep,reset);
320: if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
321: PetscCall(NEPReset_Problem(nep));
322: PetscCall(BVDestroy(&nep->V));
323: PetscCall(BVDestroy(&nep->W));
324: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
325: PetscCall(MatDestroy(&nep->resolvent));
326: nep->nwork = 0;
327: nep->state = NEP_STATE_INITIAL;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@C
332: NEPDestroy - Destroys the NEP context.
334: Collective
336: Input Parameter:
337: . nep - eigensolver context obtained from NEPCreate()
339: Level: beginner
341: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
342: @*/
343: PetscErrorCode NEPDestroy(NEP *nep)
344: {
345: PetscFunctionBegin;
346: if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
348: if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
349: PetscCall(NEPReset(*nep));
350: PetscTryTypeMethod(*nep,destroy);
351: if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
352: PetscCall(RGDestroy(&(*nep)->rg));
353: PetscCall(DSDestroy(&(*nep)->ds));
354: PetscCall(KSPDestroy(&(*nep)->refineksp));
355: PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
356: PetscCall(PetscFree((*nep)->sc));
357: /* just in case the initial vectors have not been used */
358: PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
359: if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)((*nep)->convergedctx));
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_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
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)
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: Calling sequence of fun:
663: $ PetscErrorCode fun(NEP nep,PetscScalar lambda,Mat T,Mat P,void *ctx)
664: + nep - the NEP context
665: . lambda - the scalar argument where T(.) must be evaluated
666: . T - matrix that will contain T(lambda)
667: . P - (optional) different matrix to build the preconditioner
668: - ctx - (optional) user-defined context, as set by NEPSetFunction()
670: Level: beginner
672: .seealso: NEPGetFunction(), NEPSetJacobian()
673: @*/
674: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,PetscErrorCode (*fun)(NEP nep,PetscScalar lambda,Mat T,Mat P,void *ctx),void *ctx)
675: {
676: PetscFunctionBegin;
680: if (A) PetscCheckSameComm(nep,1,A,2);
681: if (B) PetscCheckSameComm(nep,1,B,3);
683: if (nep->state) PetscCall(NEPReset(nep));
684: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
686: if (fun) nep->computefunction = fun;
687: if (ctx) nep->functionctx = ctx;
688: if (A) {
689: PetscCall(PetscObjectReference((PetscObject)A));
690: PetscCall(MatDestroy(&nep->function));
691: nep->function = A;
692: }
693: if (B) {
694: PetscCall(PetscObjectReference((PetscObject)B));
695: PetscCall(MatDestroy(&nep->function_pre));
696: nep->function_pre = B;
697: }
698: nep->fui = NEP_USER_INTERFACE_CALLBACK;
699: nep->state = NEP_STATE_INITIAL;
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: /*@C
704: NEPGetFunction - Returns the Function matrix and optionally the user
705: provided context for evaluating the Function.
707: Not Collective
709: Input Parameter:
710: . nep - the nonlinear eigensolver context
712: Output Parameters:
713: + A - location to stash Function matrix (or NULL)
714: . B - location to stash preconditioner matrix (or NULL)
715: . fun - location to put Function function (or NULL)
716: - ctx - location to stash Function context (or NULL)
718: Level: advanced
720: .seealso: NEPSetFunction()
721: @*/
722: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,PetscErrorCode (**fun)(NEP,PetscScalar,Mat,Mat,void*),void **ctx)
723: {
724: PetscFunctionBegin;
726: NEPCheckCallback(nep,1);
727: if (A) *A = nep->function;
728: if (B) *B = nep->function_pre;
729: if (fun) *fun = nep->computefunction;
730: if (ctx) *ctx = nep->functionctx;
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@C
735: NEPSetJacobian - Sets the function to compute the Jacobian T'(lambda) as well
736: as the location to store the matrix.
738: Collective
740: Input Parameters:
741: + nep - the NEP context
742: . A - Jacobian matrix
743: . jac - Jacobian evaluation routine (if NULL then NEP retains any
744: previously set value)
745: - ctx - [optional] user-defined context for private data for the Jacobian
746: evaluation routine (may be NULL) (if NULL then NEP retains any
747: previously set value)
749: Calling sequence of jac:
750: $ PetscErrorCode jac(NEP nep,PetscScalar lambda,Mat J,void *ctx)
751: + nep - the NEP context
752: . lambda - the scalar argument where T'(.) must be evaluated
753: . J - matrix that will contain T'(lambda)
754: - ctx - (optional) user-defined context, as set by NEPSetJacobian()
756: Level: beginner
758: .seealso: NEPSetFunction(), NEPGetJacobian()
759: @*/
760: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,PetscErrorCode (*jac)(NEP nep,PetscScalar lambda,Mat J,void *ctx),void *ctx)
761: {
762: PetscFunctionBegin;
765: if (A) PetscCheckSameComm(nep,1,A,2);
767: if (nep->state) PetscCall(NEPReset(nep));
768: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));
770: if (jac) nep->computejacobian = jac;
771: if (ctx) nep->jacobianctx = ctx;
772: if (A) {
773: PetscCall(PetscObjectReference((PetscObject)A));
774: PetscCall(MatDestroy(&nep->jacobian));
775: nep->jacobian = A;
776: }
777: nep->fui = NEP_USER_INTERFACE_CALLBACK;
778: nep->state = NEP_STATE_INITIAL;
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: /*@C
783: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
784: provided routine and context for evaluating the Jacobian.
786: Not Collective
788: Input Parameter:
789: . nep - the nonlinear eigensolver context
791: Output Parameters:
792: + A - location to stash Jacobian matrix (or NULL)
793: . jac - location to put Jacobian function (or NULL)
794: - ctx - location to stash Jacobian context (or NULL)
796: Level: advanced
798: .seealso: NEPSetJacobian()
799: @*/
800: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,PetscErrorCode (**jac)(NEP,PetscScalar,Mat,void*),void **ctx)
801: {
802: PetscFunctionBegin;
804: NEPCheckCallback(nep,1);
805: if (A) *A = nep->jacobian;
806: if (jac) *jac = nep->computejacobian;
807: if (ctx) *ctx = nep->jacobianctx;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: /*@
812: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
813: in split form.
815: Collective
817: Input Parameters:
818: + nep - the nonlinear eigensolver context
819: . nt - number of terms in the split form
820: . A - array of matrices
821: . f - array of functions
822: - str - structure flag for matrices
824: Notes:
825: The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
826: for i=1,...,n. The derivative T'(lambda) can be obtained using the
827: derivatives of f_i.
829: The structure flag provides information about A_i's nonzero pattern
830: (see MatStructure enum). If all matrices have the same pattern, then
831: use SAME_NONZERO_PATTERN. If the patterns are different but contained
832: in the pattern of the first one, then use SUBSET_NONZERO_PATTERN. If
833: patterns are known to be different, use DIFFERENT_NONZERO_PATTERN.
834: If set to UNKNOWN_NONZERO_PATTERN, the patterns will be compared to
835: determine if they are equal.
837: This function must be called before NEPSetUp(). If it is called again
838: after NEPSetUp() then the NEP object is reset.
840: Level: beginner
842: .seealso: NEPGetSplitOperatorTerm(), NEPGetSplitOperatorInfo(), NEPSetSplitPreconditioner()
843: @*/
844: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
845: {
846: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
848: PetscFunctionBegin;
851: PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
852: PetscAssertPointer(A,3);
853: PetscAssertPointer(f,4);
856: for (i=0;i<nt;i++) {
858: PetscCheckSameComm(nep,1,A[i],3);
860: PetscCheckSameComm(nep,1,f[i],4);
861: PetscCall(MatGetSize(A[i],&m,&n));
862: PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
863: 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);
864: 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);
865: if (!i) { m0 = m; mloc0 = mloc; }
866: 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);
867: 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);
868: PetscCall(PetscObjectReference((PetscObject)A[i]));
869: PetscCall(PetscObjectReference((PetscObject)f[i]));
870: }
872: if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
873: else PetscCall(NEPReset_Problem(nep));
875: /* allocate space and copy matrices and functions */
876: PetscCall(PetscMalloc1(nt,&nep->A));
877: for (i=0;i<nt;i++) nep->A[i] = A[i];
878: PetscCall(PetscMalloc1(nt,&nep->f));
879: for (i=0;i<nt;i++) nep->f[i] = f[i];
880: PetscCall(PetscCalloc1(nt,&nep->nrma));
881: nep->nt = nt;
882: nep->mstr = str;
883: nep->fui = NEP_USER_INTERFACE_SPLIT;
884: nep->state = NEP_STATE_INITIAL;
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
888: /*@
889: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
890: the nonlinear operator in split form.
892: Not Collective
894: Input Parameters:
895: + nep - the nonlinear eigensolver context
896: - k - the index of the requested term (starting in 0)
898: Output Parameters:
899: + A - the matrix of the requested term
900: - f - the function of the requested term
902: Level: intermediate
904: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
905: @*/
906: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
907: {
908: PetscFunctionBegin;
911: NEPCheckSplit(nep,1);
912: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
913: if (A) *A = nep->A[k];
914: if (f) *f = nep->f[k];
915: PetscFunctionReturn(PETSC_SUCCESS);
916: }
918: /*@
919: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
920: the nonlinear operator, as well as the structure flag for matrices.
922: Not Collective
924: Input Parameter:
925: . nep - the nonlinear eigensolver context
927: Output Parameters:
928: + n - the number of terms passed in NEPSetSplitOperator()
929: - str - the matrix structure flag passed in NEPSetSplitOperator()
931: Level: intermediate
933: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
934: @*/
935: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
936: {
937: PetscFunctionBegin;
939: NEPCheckSplit(nep,1);
940: if (n) *n = nep->nt;
941: if (str) *str = nep->mstr;
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: NEPSetSplitPreconditioner - Sets an operator in split form from which
947: to build the preconditioner to be used when solving the nonlinear
948: eigenvalue problem in split form.
950: Collective
952: Input Parameters:
953: + nep - the nonlinear eigensolver context
954: . ntp - number of terms in the split preconditioner
955: . P - array of matrices
956: - strp - structure flag for matrices
958: Notes:
959: The matrix for the preconditioner is expressed as P(lambda) =
960: sum_i P_i*f_i(lambda), for i=1,...,n, where the f_i functions
961: are the same as in NEPSetSplitOperator(). It is not necessary to call
962: this function. If it is not invoked, then the preconditioner is
963: built from T(lambda), i.e., both matrices and functions passed in
964: NEPSetSplitOperator().
966: The structure flag provides information about P_i's nonzero pattern
967: in the same way as in NEPSetSplitOperator().
969: If the functions defining the preconditioner operator were different
970: from the ones given in NEPSetSplitOperator(), then the split form
971: cannot be used. Use the callback interface instead.
973: Use ntp=0 to reset a previously set split preconditioner.
975: Level: advanced
977: .seealso: NEPGetSplitPreconditionerTerm(), NEPGetSplitPreconditionerInfo(), NEPSetSplitOperator()
978: @*/
979: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
980: {
981: PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
983: PetscFunctionBegin;
986: PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
987: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
988: PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
989: if (ntp) PetscAssertPointer(P,3);
992: for (i=0;i<ntp;i++) {
994: PetscCheckSameComm(nep,1,P[i],3);
995: PetscCall(MatGetSize(P[i],&m,&n));
996: PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
997: 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);
998: 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);
999: if (!i) { m0 = m; mloc0 = mloc; }
1000: 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);
1001: 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);
1002: PetscCall(PetscObjectReference((PetscObject)P[i]));
1003: }
1005: PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
1006: if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
1008: /* allocate space and copy matrices */
1009: if (ntp) {
1010: PetscCall(PetscMalloc1(ntp,&nep->P));
1011: for (i=0;i<ntp;i++) nep->P[i] = P[i];
1012: }
1013: nep->mstrp = strp;
1014: nep->state = NEP_STATE_INITIAL;
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@
1019: NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1020: the split preconditioner.
1022: Not Collective
1024: Input Parameters:
1025: + nep - the nonlinear eigensolver context
1026: - k - the index of the requested term (starting in 0)
1028: Output Parameter:
1029: . P - the matrix of the requested term
1031: Level: advanced
1033: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerInfo()
1034: @*/
1035: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1036: {
1037: PetscFunctionBegin;
1040: PetscAssertPointer(P,3);
1041: NEPCheckSplit(nep,1);
1042: PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1043: PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1044: *P = nep->P[k];
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: /*@
1049: NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1050: preconditioner, as well as the structure flag for matrices.
1052: Not Collective
1054: Input Parameter:
1055: . nep - the nonlinear eigensolver context
1057: Output Parameters:
1058: + n - the number of terms passed in NEPSetSplitPreconditioner()
1059: - strp - the matrix structure flag passed in NEPSetSplitPreconditioner()
1061: Level: advanced
1063: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerTerm()
1064: @*/
1065: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1066: {
1067: PetscFunctionBegin;
1069: NEPCheckSplit(nep,1);
1070: if (n) *n = nep->P? nep->nt: 0;
1071: if (strp) *strp = nep->mstrp;
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }