1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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> /*I "slepcnep.h" I*/
16: PetscFunctionList NEPList = 0;
17: PetscBool NEPRegisterAllCalled = PETSC_FALSE;
18: PetscClassId NEP_CLASSID = 0;
19: PetscLogEvent NEP_SetUp = 0,NEP_Solve = 0,NEP_Refine = 0,NEP_FunctionEval = 0,NEP_JacobianEval = 0,NEP_DerivativesEval = 0;
21: /*@
22: NEPCreate - Creates the default NEP context.
24: Collective on MPI_Comm
26: Input Parameter:
27: . comm - MPI communicator
29: Output Parameter:
30: . nep - location to put the NEP context
32: Level: beginner
34: .seealso: NEPSetUp(), NEPSolve(), NEPDestroy(), NEP 35: @*/
36: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep) 37: {
39: NEP nep;
43: *outnep = 0;
44: NEPInitializePackage();
45: SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView);
47: nep->max_it = 0;
48: nep->nev = 1;
49: nep->ncv = 0;
50: nep->mpd = 0;
51: nep->nini = 0;
52: nep->target = 0.0;
53: nep->tol = PETSC_DEFAULT;
54: nep->conv = NEP_CONV_REL;
55: nep->stop = NEP_STOP_BASIC;
56: nep->which = (NEPWhich)0;
57: nep->problem_type = (NEPProblemType)0;
58: nep->refine = NEP_REFINE_NONE;
59: nep->npart = 1;
60: nep->rtol = PETSC_DEFAULT;
61: nep->rits = PETSC_DEFAULT;
62: nep->scheme = (NEPRefineScheme)0;
63: nep->trackall = PETSC_FALSE;
65: nep->computefunction = NULL;
66: nep->computejacobian = NULL;
67: nep->functionctx = NULL;
68: nep->jacobianctx = NULL;
69: nep->computederivatives = NULL;
70: nep->derivativesctx = NULL;
71: nep->converged = NEPConvergedRelative;
72: nep->convergeduser = NULL;
73: nep->convergeddestroy= NULL;
74: nep->stopping = NEPStoppingBasic;
75: nep->stoppinguser = NULL;
76: nep->stoppingdestroy = NULL;
77: nep->convergedctx = NULL;
78: nep->stoppingctx = NULL;
79: nep->numbermonitors = 0;
81: nep->ds = NULL;
82: nep->V = NULL;
83: nep->rg = NULL;
84: nep->function = NULL;
85: nep->function_pre = NULL;
86: nep->jacobian = NULL;
87: nep->derivatives = NULL;
88: nep->A = NULL;
89: nep->f = NULL;
90: nep->nt = 0;
91: nep->mstr = DIFFERENT_NONZERO_PATTERN;
92: nep->IS = NULL;
93: nep->eigr = NULL;
94: nep->eigi = NULL;
95: nep->errest = NULL;
96: nep->perm = NULL;
97: nep->nwork = 0;
98: nep->work = NULL;
99: nep->data = NULL;
101: nep->state = NEP_STATE_INITIAL;
102: nep->nconv = 0;
103: nep->its = 0;
104: nep->n = 0;
105: nep->nloc = 0;
106: nep->nrma = NULL;
107: nep->fui = (NEPUserInterface)0;
108: nep->reason = NEP_CONVERGED_ITERATING;
110: PetscNewLog(nep,&nep->sc);
111: *outnep = nep;
112: return(0);
113: }
115: /*@C
116: NEPSetType - Selects the particular solver to be used in the NEP object.
118: Logically Collective on NEP120: Input Parameters:
121: + nep - the nonlinear eigensolver context
122: - type - a known method
124: Options Database Key:
125: . -nep_type <method> - Sets the method; use -help for a list
126: of available methods
128: Notes:
129: See "slepc/include/slepcnep.h" for available methods.
131: Normally, it is best to use the NEPSetFromOptions() command and
132: then set the NEP type from the options database rather than by using
133: this routine. Using the options database provides the user with
134: maximum flexibility in evaluating the different available methods.
135: The NEPSetType() routine is provided for those situations where it
136: is necessary to set the iterative solver independently of the command
137: line or options database.
139: Level: intermediate
141: .seealso: NEPType142: @*/
143: PetscErrorCode NEPSetType(NEP nep,NEPType type)144: {
145: PetscErrorCode ierr,(*r)(NEP);
146: PetscBool match;
152: PetscObjectTypeCompare((PetscObject)nep,type,&match);
153: if (match) return(0);
155: PetscFunctionListFind(NEPList,type,&r);
156: if (!r) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);
158: if (nep->ops->destroy) { (*nep->ops->destroy)(nep); }
159: PetscMemzero(nep->ops,sizeof(struct _NEPOps));
161: nep->state = NEP_STATE_INITIAL;
162: PetscObjectChangeTypeName((PetscObject)nep,type);
163: (*r)(nep);
164: return(0);
165: }
167: /*@C
168: NEPGetType - Gets the NEP type as a string from the NEP object.
170: Not Collective
172: Input Parameter:
173: . nep - the eigensolver context
175: Output Parameter:
176: . name - name of NEP method
178: Level: intermediate
180: .seealso: NEPSetType()
181: @*/
182: PetscErrorCode NEPGetType(NEP nep,NEPType *type)183: {
187: *type = ((PetscObject)nep)->type_name;
188: return(0);
189: }
191: /*@C
192: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
194: Not Collective
196: Input Parameters:
197: + name - name of a new user-defined solver
198: - function - routine to create the solver context
200: Notes:
201: NEPRegister() may be called multiple times to add several user-defined solvers.
203: Sample usage:
204: .vb
205: NEPRegister("my_solver",MySolverCreate);
206: .ve
208: Then, your solver can be chosen with the procedural interface via
209: $ NEPSetType(nep,"my_solver")
210: or at runtime via the option
211: $ -nep_type my_solver
213: Level: advanced
215: .seealso: NEPRegisterAll()
216: @*/
217: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))218: {
222: NEPInitializePackage();
223: PetscFunctionListAdd(&NEPList,name,function);
224: return(0);
225: }
227: /*
228: NEPReset_Problem - Destroys the problem matrices.
229: @*/
230: PetscErrorCode NEPReset_Problem(NEP nep)231: {
233: PetscInt i;
237: MatDestroy(&nep->function);
238: MatDestroy(&nep->function_pre);
239: MatDestroy(&nep->jacobian);
240: MatDestroy(&nep->derivatives);
241: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
242: MatDestroyMatrices(nep->nt,&nep->A);
243: for (i=0;i<nep->nt;i++) {
244: FNDestroy(&nep->f[i]);
245: }
246: PetscFree(nep->f);
247: PetscFree(nep->nrma);
248: nep->nt = 0;
249: }
250: return(0);
251: }
252: /*@
253: NEPReset - Resets the NEP context to the initial state (prior to setup)
254: and destroys any allocated Vecs and Mats.
256: Collective on NEP258: Input Parameter:
259: . nep - eigensolver context obtained from NEPCreate()
261: Level: advanced
263: .seealso: NEPDestroy()
264: @*/
265: PetscErrorCode NEPReset(NEP nep)266: {
271: if (!nep) return(0);
272: if (nep->ops->reset) { (nep->ops->reset)(nep); }
273: if (nep->refineksp) { KSPReset(nep->refineksp); }
274: NEPReset_Problem(nep);
275: BVDestroy(&nep->V);
276: VecDestroyVecs(nep->nwork,&nep->work);
277: nep->nwork = 0;
278: nep->state = NEP_STATE_INITIAL;
279: return(0);
280: }
282: /*@
283: NEPDestroy - Destroys the NEP context.
285: Collective on NEP287: Input Parameter:
288: . nep - eigensolver context obtained from NEPCreate()
290: Level: beginner
292: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
293: @*/
294: PetscErrorCode NEPDestroy(NEP *nep)295: {
299: if (!*nep) return(0);
301: if (--((PetscObject)(*nep))->refct > 0) { *nep = 0; return(0); }
302: NEPReset(*nep);
303: if ((*nep)->ops->destroy) { (*(*nep)->ops->destroy)(*nep); }
304: if ((*nep)->eigr) {
305: PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm);
306: }
307: RGDestroy(&(*nep)->rg);
308: DSDestroy(&(*nep)->ds);
309: KSPDestroy(&(*nep)->refineksp);
310: PetscSubcommDestroy(&(*nep)->refinesubc);
311: PetscFree((*nep)->sc);
312: /* just in case the initial vectors have not been used */
313: SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS);
314: if ((*nep)->convergeddestroy) {
315: (*(*nep)->convergeddestroy)((*nep)->convergedctx);
316: }
317: NEPMonitorCancel(*nep);
318: PetscHeaderDestroy(nep);
319: return(0);
320: }
322: /*@
323: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
325: Collective on NEP327: Input Parameters:
328: + nep - eigensolver context obtained from NEPCreate()
329: - bv - the basis vectors object
331: Note:
332: Use NEPGetBV() to retrieve the basis vectors context (for example,
333: to free it at the end of the computations).
335: Level: advanced
337: .seealso: NEPGetBV()
338: @*/
339: PetscErrorCode NEPSetBV(NEP nep,BV bv)340: {
347: PetscObjectReference((PetscObject)bv);
348: BVDestroy(&nep->V);
349: nep->V = bv;
350: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
351: return(0);
352: }
354: /*@
355: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
356: eigensolver object.
358: Not Collective
360: Input Parameters:
361: . nep - eigensolver context obtained from NEPCreate()
363: Output Parameter:
364: . bv - basis vectors context
366: Level: advanced
368: .seealso: NEPSetBV()
369: @*/
370: PetscErrorCode NEPGetBV(NEP nep,BV *bv)371: {
377: if (!nep->V) {
378: BVCreate(PetscObjectComm((PetscObject)nep),&nep->V);
379: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
380: }
381: *bv = nep->V;
382: return(0);
383: }
385: /*@
386: NEPSetRG - Associates a region object to the nonlinear eigensolver.
388: Collective on NEP390: Input Parameters:
391: + nep - eigensolver context obtained from NEPCreate()
392: - rg - the region object
394: Note:
395: Use NEPGetRG() to retrieve the region context (for example,
396: to free it at the end of the computations).
398: Level: advanced
400: .seealso: NEPGetRG()
401: @*/
402: PetscErrorCode NEPSetRG(NEP nep,RG rg)403: {
410: PetscObjectReference((PetscObject)rg);
411: RGDestroy(&nep->rg);
412: nep->rg = rg;
413: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
414: return(0);
415: }
417: /*@
418: NEPGetRG - Obtain the region object associated to the
419: nonlinear eigensolver object.
421: Not Collective
423: Input Parameters:
424: . nep - eigensolver context obtained from NEPCreate()
426: Output Parameter:
427: . rg - region context
429: Level: advanced
431: .seealso: NEPSetRG()
432: @*/
433: PetscErrorCode NEPGetRG(NEP nep,RG *rg)434: {
440: if (!nep->rg) {
441: RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg);
442: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
443: }
444: *rg = nep->rg;
445: return(0);
446: }
448: /*@
449: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
451: Collective on NEP453: Input Parameters:
454: + nep - eigensolver context obtained from NEPCreate()
455: - ds - the direct solver object
457: Note:
458: Use NEPGetDS() to retrieve the direct solver context (for example,
459: to free it at the end of the computations).
461: Level: advanced
463: .seealso: NEPGetDS()
464: @*/
465: PetscErrorCode NEPSetDS(NEP nep,DS ds)466: {
473: PetscObjectReference((PetscObject)ds);
474: DSDestroy(&nep->ds);
475: nep->ds = ds;
476: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
477: return(0);
478: }
480: /*@
481: NEPGetDS - Obtain the direct solver object associated to the
482: nonlinear eigensolver object.
484: Not Collective
486: Input Parameters:
487: . nep - eigensolver context obtained from NEPCreate()
489: Output Parameter:
490: . ds - direct solver context
492: Level: advanced
494: .seealso: NEPSetDS()
495: @*/
496: PetscErrorCode NEPGetDS(NEP nep,DS *ds)497: {
503: if (!nep->ds) {
504: DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds);
505: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
506: }
507: *ds = nep->ds;
508: return(0);
509: }
511: /*@
512: NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
513: object in the refinement phase.
515: Not Collective
517: Input Parameters:
518: . nep - eigensolver context obtained from NEPCreate()
520: Output Parameter:
521: . ksp - ksp context
523: Level: advanced
525: .seealso: NEPSetRefine()
526: @*/
527: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)528: {
534: if (!nep->refineksp) {
535: if (nep->npart>1) {
536: /* Split in subcomunicators */
537: PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc);
538: PetscSubcommSetNumber(nep->refinesubc,nep->npart);
539: PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
540: PetscLogObjectMemory((PetscObject)nep,sizeof(PetscSubcomm));
541: }
542: KSPCreate((nep->npart==1)?PetscObjectComm((PetscObject)nep):PetscSubcommChild(nep->refinesubc),&nep->refineksp);
543: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->refineksp);
544: KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix);
545: KSPAppendOptionsPrefix(*ksp,"nep_refine_");
546: KSPSetTolerances(nep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
547: }
548: *ksp = nep->refineksp;
549: return(0);
550: }
552: /*@
553: NEPSetTarget - Sets the value of the target.
555: Logically Collective on NEP557: Input Parameters:
558: + nep - eigensolver context
559: - target - the value of the target
561: Options Database Key:
562: . -nep_target <scalar> - the value of the target
564: Notes:
565: The target is a scalar value used to determine the portion of the spectrum
566: of interest. It is used in combination with NEPSetWhichEigenpairs().
568: In the case of complex scalars, a complex value can be provided in the
569: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
570: -nep_target 1.0+2.0i
572: Level: intermediate
574: .seealso: NEPGetTarget(), NEPSetWhichEigenpairs()
575: @*/
576: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)577: {
581: nep->target = target;
582: return(0);
583: }
585: /*@
586: NEPGetTarget - Gets the value of the target.
588: Not Collective
590: Input Parameter:
591: . nep - eigensolver context
593: Output Parameter:
594: . target - the value of the target
596: Note:
597: If the target was not set by the user, then zero is returned.
599: Level: intermediate
601: .seealso: NEPSetTarget()
602: @*/
603: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)604: {
608: *target = nep->target;
609: return(0);
610: }
612: /*@C
613: NEPSetFunction - Sets the function to compute the nonlinear Function T(lambda)
614: as well as the location to store the matrix.
616: Logically Collective on NEP and Mat
618: Input Parameters:
619: + nep - the NEP context
620: . A - Function matrix
621: . B - preconditioner matrix (usually same as the Function)
622: . fun - Function evaluation routine (if NULL then NEP retains any
623: previously set value)
624: - ctx - [optional] user-defined context for private data for the Function
625: evaluation routine (may be NULL) (if NULL then NEP retains any
626: previously set value)
628: Calling Sequence of fun:
629: $ fun(NEP nep,PetscScalar lambda,Mat F,Mat P,void *ctx)
631: + nep - the NEP context
632: . lambda - the scalar argument where T(.) must be evaluated
633: . T - matrix that will contain T(lambda)
634: . P - (optional) different matrix to build the preconditioner
635: - ctx - (optional) user-defined context, as set by NEPSetFunction()
637: Level: beginner
639: .seealso: NEPGetFunction(), NEPSetJacobian()
640: @*/
641: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,PetscErrorCode (*fun)(NEP,PetscScalar,Mat,Mat,void*),void *ctx)642: {
652: if (nep->state) { NEPReset(nep); }
653: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }
655: if (fun) nep->computefunction = fun;
656: if (ctx) nep->functionctx = ctx;
657: if (A) {
658: PetscObjectReference((PetscObject)A);
659: MatDestroy(&nep->function);
660: nep->function = A;
661: }
662: if (B) {
663: PetscObjectReference((PetscObject)B);
664: MatDestroy(&nep->function_pre);
665: nep->function_pre = B;
666: }
667: nep->fui = NEP_USER_INTERFACE_CALLBACK;
668: nep->state = NEP_STATE_INITIAL;
669: return(0);
670: }
672: /*@C
673: NEPGetFunction - Returns the Function matrix and optionally the user
674: provided context for evaluating the Function.
676: Not Collective, but Mat object will be parallel if NEP object is
678: Input Parameter:
679: . nep - the nonlinear eigensolver context
681: Output Parameters:
682: + A - location to stash Function matrix (or NULL)
683: . B - location to stash preconditioner matrix (or NULL)
684: . fun - location to put Function function (or NULL)
685: - ctx - location to stash Function context (or NULL)
687: Level: advanced
689: .seealso: NEPSetFunction()
690: @*/
691: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,PetscErrorCode (**fun)(NEP,PetscScalar,Mat,Mat,void*),void **ctx)692: {
695: NEPCheckCallback(nep,1);
696: if (A) *A = nep->function;
697: if (B) *B = nep->function_pre;
698: if (fun) *fun = nep->computefunction;
699: if (ctx) *ctx = nep->functionctx;
700: return(0);
701: }
703: /*@C
704: NEPSetJacobian - Sets the function to compute Jacobian T'(lambda) as well
705: as the location to store the matrix.
707: Logically Collective on NEP and Mat
709: Input Parameters:
710: + nep - the NEP context
711: . A - Jacobian matrix
712: . jac - Jacobian evaluation routine (if NULL then NEP retains any
713: previously set value)
714: - ctx - [optional] user-defined context for private data for the Jacobian
715: evaluation routine (may be NULL) (if NULL then NEP retains any
716: previously set value)
718: Calling Sequence of jac:
719: $ jac(NEP nep,PetscScalar lambda,Mat J,void *ctx)
721: + nep - the NEP context
722: . lambda - the scalar argument where T'(.) must be evaluated
723: . J - matrix that will contain T'(lambda)
724: - ctx - (optional) user-defined context, as set by NEPSetJacobian()
726: Level: beginner
728: .seealso: NEPSetFunction(), NEPGetJacobian()
729: @*/
730: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,PetscErrorCode (*jac)(NEP,PetscScalar,Mat,void*),void *ctx)731: {
739: if (nep->state) { NEPReset(nep); }
740: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }
742: if (jac) nep->computejacobian = jac;
743: if (ctx) nep->jacobianctx = ctx;
744: if (A) {
745: PetscObjectReference((PetscObject)A);
746: MatDestroy(&nep->jacobian);
747: nep->jacobian = A;
748: }
749: nep->fui = NEP_USER_INTERFACE_CALLBACK;
750: nep->state = NEP_STATE_INITIAL;
751: return(0);
752: }
754: /*@C
755: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
756: provided routine and context for evaluating the Jacobian.
758: Not Collective, but Mat object will be parallel if NEP object is
760: Input Parameter:
761: . nep - the nonlinear eigensolver context
763: Output Parameters:
764: + A - location to stash Jacobian matrix (or NULL)
765: . jac - location to put Jacobian function (or NULL)
766: - ctx - location to stash Jacobian context (or NULL)
768: Level: advanced
770: .seealso: NEPSetJacobian()
771: @*/
772: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,PetscErrorCode (**jac)(NEP,PetscScalar,Mat,void*),void **ctx)773: {
776: NEPCheckCallback(nep,1);
777: if (A) *A = nep->jacobian;
778: if (jac) *jac = nep->computejacobian;
779: if (ctx) *ctx = nep->jacobianctx;
780: return(0);
781: }
783: /*@
784: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
785: in split form.
787: Collective on NEP, Mat and FN789: Input Parameters:
790: + nep - the nonlinear eigensolver context
791: . n - number of terms in the split form
792: . A - array of matrices
793: . f - array of functions
794: - str - structure flag for matrices
796: Notes:
797: The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
798: for i=1,...,n. The derivative T'(lambda) can be obtained using the
799: derivatives of f_i.
801: The structure flag provides information about A_i's nonzero pattern
802: (see MatStructure enum). If all matrices have the same pattern, then
803: use SAME_NONZERO_PATTERN. If the patterns are different but contained
804: in the pattern of the first one, then use SUBSET_NONZERO_PATTERN.
805: Otherwise use DIFFERENT_NONZERO_PATTERN.
807: This function must be called before NEPSetUp(). If it is called again
808: after NEPSetUp() then the NEP object is reset.
810: Level: beginner
812: .seealso: NEPGetSplitOperatorTerm(), NEPGetSplitOperatorInfo()
813: @*/
814: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt n,Mat A[],FN f[],MatStructure str)815: {
816: PetscInt i;
822: if (n <= 0) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %D",n);
828: for (i=0;i<n;i++) {
831: PetscObjectReference((PetscObject)A[i]);
832: PetscObjectReference((PetscObject)f[i]);
833: }
835: if (nep->state) { NEPReset(nep); }
836: else { NEPReset_Problem(nep); }
838: /* allocate space and copy matrices and functions */
839: PetscMalloc1(n,&nep->A);
840: PetscLogObjectMemory((PetscObject)nep,n*sizeof(Mat));
841: for (i=0;i<n;i++) nep->A[i] = A[i];
842: PetscMalloc1(n,&nep->f);
843: PetscLogObjectMemory((PetscObject)nep,n*sizeof(FN));
844: for (i=0;i<n;i++) nep->f[i] = f[i];
845: PetscCalloc1(n,&nep->nrma);
846: PetscLogObjectMemory((PetscObject)nep,n*sizeof(PetscReal));
847: nep->nt = n;
848: nep->mstr = str;
849: nep->fui = NEP_USER_INTERFACE_SPLIT;
850: nep->state = NEP_STATE_INITIAL;
851: return(0);
852: }
854: /*@
855: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
856: the nonlinear operator in split form.
858: Not collective, though parallel Mats and FNs are returned if the NEP is parallel
860: Input Parameter:
861: + nep - the nonlinear eigensolver context
862: - k - the index of the requested term (starting in 0)
864: Output Parameters:
865: + A - the matrix of the requested term
866: - f - the function of the requested term
868: Level: intermediate
870: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
871: @*/
872: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)873: {
876: NEPCheckSplit(nep,1);
877: if (k<0 || k>=nep->nt) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",nep->nt-1);
878: if (A) *A = nep->A[k];
879: if (f) *f = nep->f[k];
880: return(0);
881: }
883: /*@
884: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
885: the nonlinear operator, as well as the structure flag for matrices.
887: Not collective
889: Input Parameter:
890: . nep - the nonlinear eigensolver context
892: Output Parameters:
893: + n - the number of terms passed in NEPSetSplitOperator()
894: - str - the matrix structure flag passed in NEPSetSplitOperator()
896: Level: intermediate
898: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
899: @*/
900: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)901: {
904: NEPCheckSplit(nep,1);
905: if (n) *n = nep->nt;
906: if (str) *str = nep->mstr;
907: return(0);
908: }
910: /*@C
911: NEPSetDerivatives - Sets the function to compute the k-th derivative T^(k)(lambda)
912: for any value of k (including 0), as well as the location to store the matrix.
914: Logically Collective on NEP and Mat
916: Input Parameters:
917: + nep - the NEP context
918: . A - the matrix to store the computed derivative
919: . der - routing to evaluate the k-th derivative (if NULL then NEP retains any
920: previously set value)
921: - ctx - [optional] user-defined context for private data for the derivatives
922: evaluation routine (may be NULL) (if NULL then NEP retains any
923: previously set value)
925: Level: beginner
927: .seealso: NEPSetFunction(), NEPGetDerivatives()
928: @*/
929: PetscErrorCode NEPSetDerivatives(NEP nep,Mat A,PetscErrorCode (*der)(NEP,PetscScalar,PetscInt,Mat,void*),void *ctx)930: {
938: if (nep->state) { NEPReset(nep); }
939: else { NEPReset_Problem(nep); }
941: if (der) nep->computederivatives = der;
942: if (ctx) nep->derivativesctx = ctx;
943: if (A) {
944: PetscObjectReference((PetscObject)A);
945: MatDestroy(&nep->derivatives);
946: nep->derivatives = A;
947: }
948: nep->fui = NEP_USER_INTERFACE_DERIVATIVES;
949: nep->state = NEP_STATE_INITIAL;
950: return(0);
951: }
953: /*@C
954: NEPGetDerivatives - Returns the derivatives matrix and optionally the user
955: provided routine and context for evaluating the derivatives.
957: Not Collective, but Mat object will be parallel if NEP object is
959: Input Parameter:
960: . nep - the nonlinear eigensolver context
962: Output Parameters:
963: + A - location to stash the derivatives matrix (or NULL)
964: . der - location to put derivatives function (or NULL)
965: - ctx - location to stash derivatives context (or NULL)
967: Level: advanced
969: .seealso: NEPSetDerivatives()
970: @*/
971: PetscErrorCode NEPGetDerivatives(NEP nep,Mat *A,PetscErrorCode (**der)(NEP,PetscScalar,PetscInt,Mat,void*),void **ctx)972: {
975: NEPCheckDerivatives(nep,1);
976: if (A) *A = nep->derivatives;
977: if (der) *der = nep->computederivatives;
978: if (ctx) *ctx = nep->derivativesctx;
979: return(0);
980: }