Actual source code: nepbasic.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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,NEP_Resolvent = 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;
64: nep->twosided = PETSC_FALSE;
66: nep->computefunction = NULL;
67: nep->computejacobian = NULL;
68: nep->functionctx = NULL;
69: nep->jacobianctx = NULL;
70: nep->computederivatives = NULL;
71: nep->derivativesctx = NULL;
72: nep->converged = NEPConvergedRelative;
73: nep->convergeduser = NULL;
74: nep->convergeddestroy= NULL;
75: nep->stopping = NEPStoppingBasic;
76: nep->stoppinguser = NULL;
77: nep->stoppingdestroy = NULL;
78: nep->convergedctx = NULL;
79: nep->stoppingctx = NULL;
80: nep->numbermonitors = 0;
82: nep->ds = NULL;
83: nep->V = NULL;
84: nep->W = NULL;
85: nep->rg = NULL;
86: nep->function = NULL;
87: nep->function_pre = NULL;
88: nep->jacobian = NULL;
89: nep->derivatives = NULL;
90: nep->A = NULL;
91: nep->f = NULL;
92: nep->nt = 0;
93: nep->mstr = DIFFERENT_NONZERO_PATTERN;
94: nep->IS = NULL;
95: nep->eigr = NULL;
96: nep->eigi = NULL;
97: nep->errest = NULL;
98: nep->perm = NULL;
99: nep->nwork = 0;
100: nep->work = NULL;
101: nep->data = NULL;
103: nep->state = NEP_STATE_INITIAL;
104: nep->nconv = 0;
105: nep->its = 0;
106: nep->n = 0;
107: nep->nloc = 0;
108: nep->nrma = NULL;
109: nep->fui = (NEPUserInterface)0;
110: nep->useds = PETSC_FALSE;
111: nep->hasts = PETSC_FALSE;
112: nep->resolvent = NULL;
113: nep->reason = NEP_CONVERGED_ITERATING;
115: PetscNewLog(nep,&nep->sc);
116: *outnep = nep;
117: return(0);
118: }
120: /*@C
121: NEPSetType - Selects the particular solver to be used in the NEP object.
123: Logically Collective on NEP
125: Input Parameters:
126: + nep - the nonlinear eigensolver context
127: - type - a known method
129: Options Database Key:
130: . -nep_type <method> - Sets the method; use -help for a list
131: of available methods
133: Notes:
134: See "slepc/include/slepcnep.h" for available methods.
136: Normally, it is best to use the NEPSetFromOptions() command and
137: then set the NEP type from the options database rather than by using
138: this routine. Using the options database provides the user with
139: maximum flexibility in evaluating the different available methods.
140: The NEPSetType() routine is provided for those situations where it
141: is necessary to set the iterative solver independently of the command
142: line or options database.
144: Level: intermediate
146: .seealso: NEPType
147: @*/
148: PetscErrorCode NEPSetType(NEP nep,NEPType type)
149: {
150: PetscErrorCode ierr,(*r)(NEP);
151: PetscBool match;
157: PetscObjectTypeCompare((PetscObject)nep,type,&match);
158: if (match) return(0);
160: PetscFunctionListFind(NEPList,type,&r);
161: if (!r) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);
163: if (nep->ops->destroy) { (*nep->ops->destroy)(nep); }
164: PetscMemzero(nep->ops,sizeof(struct _NEPOps));
166: nep->state = NEP_STATE_INITIAL;
167: PetscObjectChangeTypeName((PetscObject)nep,type);
168: (*r)(nep);
169: return(0);
170: }
172: /*@C
173: NEPGetType - Gets the NEP type as a string from the NEP object.
175: Not Collective
177: Input Parameter:
178: . nep - the eigensolver context
180: Output Parameter:
181: . name - name of NEP method
183: Level: intermediate
185: .seealso: NEPSetType()
186: @*/
187: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
188: {
192: *type = ((PetscObject)nep)->type_name;
193: return(0);
194: }
196: /*@C
197: NEPRegister - Adds a method to the nonlinear eigenproblem solver package.
199: Not Collective
201: Input Parameters:
202: + name - name of a new user-defined solver
203: - function - routine to create the solver context
205: Notes:
206: NEPRegister() may be called multiple times to add several user-defined solvers.
208: Sample usage:
209: .vb
210: NEPRegister("my_solver",MySolverCreate);
211: .ve
213: Then, your solver can be chosen with the procedural interface via
214: $ NEPSetType(nep,"my_solver")
215: or at runtime via the option
216: $ -nep_type my_solver
218: Level: advanced
220: .seealso: NEPRegisterAll()
221: @*/
222: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
223: {
227: NEPInitializePackage();
228: PetscFunctionListAdd(&NEPList,name,function);
229: return(0);
230: }
232: /*
233: NEPReset_Problem - Destroys the problem matrices.
234: @*/
235: PetscErrorCode NEPReset_Problem(NEP nep)
236: {
238: PetscInt i;
242: MatDestroy(&nep->function);
243: MatDestroy(&nep->function_pre);
244: MatDestroy(&nep->jacobian);
245: MatDestroy(&nep->derivatives);
246: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
247: MatDestroyMatrices(nep->nt,&nep->A);
248: for (i=0;i<nep->nt;i++) {
249: FNDestroy(&nep->f[i]);
250: }
251: PetscFree(nep->f);
252: PetscFree(nep->nrma);
253: nep->nt = 0;
254: }
255: return(0);
256: }
257: /*@
258: NEPReset - Resets the NEP context to the initial state (prior to setup)
259: and destroys any allocated Vecs and Mats.
261: Collective on NEP
263: Input Parameter:
264: . nep - eigensolver context obtained from NEPCreate()
266: Level: advanced
268: .seealso: NEPDestroy()
269: @*/
270: PetscErrorCode NEPReset(NEP nep)
271: {
276: if (!nep) return(0);
277: if (nep->ops->reset) { (nep->ops->reset)(nep); }
278: if (nep->refineksp) { KSPReset(nep->refineksp); }
279: NEPReset_Problem(nep);
280: BVDestroy(&nep->V);
281: BVDestroy(&nep->W);
282: VecDestroyVecs(nep->nwork,&nep->work);
283: MatDestroy(&nep->resolvent);
284: nep->nwork = 0;
285: nep->state = NEP_STATE_INITIAL;
286: return(0);
287: }
289: /*@
290: NEPDestroy - Destroys the NEP context.
292: Collective on NEP
294: Input Parameter:
295: . nep - eigensolver context obtained from NEPCreate()
297: Level: beginner
299: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
300: @*/
301: PetscErrorCode NEPDestroy(NEP *nep)
302: {
306: if (!*nep) return(0);
308: if (--((PetscObject)(*nep))->refct > 0) { *nep = 0; return(0); }
309: NEPReset(*nep);
310: if ((*nep)->ops->destroy) { (*(*nep)->ops->destroy)(*nep); }
311: if ((*nep)->eigr) {
312: PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm);
313: }
314: RGDestroy(&(*nep)->rg);
315: DSDestroy(&(*nep)->ds);
316: KSPDestroy(&(*nep)->refineksp);
317: PetscSubcommDestroy(&(*nep)->refinesubc);
318: PetscFree((*nep)->sc);
319: /* just in case the initial vectors have not been used */
320: SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS);
321: if ((*nep)->convergeddestroy) {
322: (*(*nep)->convergeddestroy)((*nep)->convergedctx);
323: }
324: NEPMonitorCancel(*nep);
325: PetscHeaderDestroy(nep);
326: return(0);
327: }
329: /*@
330: NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.
332: Collective on NEP
334: Input Parameters:
335: + nep - eigensolver context obtained from NEPCreate()
336: - bv - the basis vectors object
338: Note:
339: Use NEPGetBV() to retrieve the basis vectors context (for example,
340: to free it at the end of the computations).
342: Level: advanced
344: .seealso: NEPGetBV()
345: @*/
346: PetscErrorCode NEPSetBV(NEP nep,BV bv)
347: {
354: PetscObjectReference((PetscObject)bv);
355: BVDestroy(&nep->V);
356: nep->V = bv;
357: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
358: return(0);
359: }
361: /*@
362: NEPGetBV - Obtain the basis vectors object associated to the nonlinear
363: eigensolver object.
365: Not Collective
367: Input Parameters:
368: . nep - eigensolver context obtained from NEPCreate()
370: Output Parameter:
371: . bv - basis vectors context
373: Level: advanced
375: .seealso: NEPSetBV()
376: @*/
377: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
378: {
384: if (!nep->V) {
385: BVCreate(PetscObjectComm((PetscObject)nep),&nep->V);
386: PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0);
387: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
388: PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options);
389: }
390: *bv = nep->V;
391: return(0);
392: }
394: /*@
395: NEPSetRG - Associates a region object to the nonlinear eigensolver.
397: Collective on NEP
399: Input Parameters:
400: + nep - eigensolver context obtained from NEPCreate()
401: - rg - the region object
403: Note:
404: Use NEPGetRG() to retrieve the region context (for example,
405: to free it at the end of the computations).
407: Level: advanced
409: .seealso: NEPGetRG()
410: @*/
411: PetscErrorCode NEPSetRG(NEP nep,RG rg)
412: {
419: PetscObjectReference((PetscObject)rg);
420: RGDestroy(&nep->rg);
421: nep->rg = rg;
422: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
423: return(0);
424: }
426: /*@
427: NEPGetRG - Obtain the region object associated to the
428: nonlinear eigensolver object.
430: Not Collective
432: Input Parameters:
433: . nep - eigensolver context obtained from NEPCreate()
435: Output Parameter:
436: . rg - region context
438: Level: advanced
440: .seealso: NEPSetRG()
441: @*/
442: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
443: {
449: if (!nep->rg) {
450: RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg);
451: PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0);
452: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
453: PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options);
454: }
455: *rg = nep->rg;
456: return(0);
457: }
459: /*@
460: NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.
462: Collective on NEP
464: Input Parameters:
465: + nep - eigensolver context obtained from NEPCreate()
466: - ds - the direct solver object
468: Note:
469: Use NEPGetDS() to retrieve the direct solver context (for example,
470: to free it at the end of the computations).
472: Level: advanced
474: .seealso: NEPGetDS()
475: @*/
476: PetscErrorCode NEPSetDS(NEP nep,DS ds)
477: {
484: PetscObjectReference((PetscObject)ds);
485: DSDestroy(&nep->ds);
486: nep->ds = ds;
487: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
488: return(0);
489: }
491: /*@
492: NEPGetDS - Obtain the direct solver object associated to the
493: nonlinear eigensolver object.
495: Not Collective
497: Input Parameters:
498: . nep - eigensolver context obtained from NEPCreate()
500: Output Parameter:
501: . ds - direct solver context
503: Level: advanced
505: .seealso: NEPSetDS()
506: @*/
507: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
508: {
514: if (!nep->ds) {
515: DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds);
516: PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0);
517: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
518: PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options);
519: }
520: *ds = nep->ds;
521: return(0);
522: }
524: /*@
525: NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
526: object in the refinement phase.
528: Not Collective
530: Input Parameters:
531: . nep - eigensolver context obtained from NEPCreate()
533: Output Parameter:
534: . ksp - ksp context
536: Level: advanced
538: .seealso: NEPSetRefine()
539: @*/
540: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
541: {
547: if (!nep->refineksp) {
548: if (nep->npart>1) {
549: /* Split in subcomunicators */
550: PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc);
551: PetscSubcommSetNumber(nep->refinesubc,nep->npart);
552: PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
553: PetscLogObjectMemory((PetscObject)nep,sizeof(PetscSubcomm));
554: }
555: KSPCreate((nep->npart==1)?PetscObjectComm((PetscObject)nep):PetscSubcommChild(nep->refinesubc),&nep->refineksp);
556: PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0);
557: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->refineksp);
558: PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options);
559: KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix);
560: KSPAppendOptionsPrefix(*ksp,"nep_refine_");
561: KSPSetTolerances(nep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
562: }
563: *ksp = nep->refineksp;
564: return(0);
565: }
567: /*@
568: NEPSetTarget - Sets the value of the target.
570: Logically Collective on NEP
572: Input Parameters:
573: + nep - eigensolver context
574: - target - the value of the target
576: Options Database Key:
577: . -nep_target <scalar> - the value of the target
579: Notes:
580: The target is a scalar value used to determine the portion of the spectrum
581: of interest. It is used in combination with NEPSetWhichEigenpairs().
583: In the case of complex scalars, a complex value can be provided in the
584: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
585: -nep_target 1.0+2.0i
587: Level: intermediate
589: .seealso: NEPGetTarget(), NEPSetWhichEigenpairs()
590: @*/
591: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
592: {
596: nep->target = target;
597: return(0);
598: }
600: /*@
601: NEPGetTarget - Gets the value of the target.
603: Not Collective
605: Input Parameter:
606: . nep - eigensolver context
608: Output Parameter:
609: . target - the value of the target
611: Note:
612: If the target was not set by the user, then zero is returned.
614: Level: intermediate
616: .seealso: NEPSetTarget()
617: @*/
618: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
619: {
623: *target = nep->target;
624: return(0);
625: }
627: /*@C
628: NEPSetFunction - Sets the function to compute the nonlinear Function T(lambda)
629: as well as the location to store the matrix.
631: Logically Collective on NEP and Mat
633: Input Parameters:
634: + nep - the NEP context
635: . A - Function matrix
636: . B - preconditioner matrix (usually same as the Function)
637: . fun - Function evaluation routine (if NULL then NEP retains any
638: previously set value)
639: - ctx - [optional] user-defined context for private data for the Function
640: evaluation routine (may be NULL) (if NULL then NEP retains any
641: previously set value)
643: Calling Sequence of fun:
644: $ fun(NEP nep,PetscScalar lambda,Mat F,Mat P,void *ctx)
646: + nep - the NEP context
647: . lambda - the scalar argument where T(.) must be evaluated
648: . T - matrix that will contain T(lambda)
649: . P - (optional) different matrix to build the preconditioner
650: - ctx - (optional) user-defined context, as set by NEPSetFunction()
652: Level: beginner
654: .seealso: NEPGetFunction(), NEPSetJacobian()
655: @*/
656: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,PetscErrorCode (*fun)(NEP,PetscScalar,Mat,Mat,void*),void *ctx)
657: {
667: if (nep->state) { NEPReset(nep); }
668: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }
670: if (fun) nep->computefunction = fun;
671: if (ctx) nep->functionctx = ctx;
672: if (A) {
673: PetscObjectReference((PetscObject)A);
674: MatDestroy(&nep->function);
675: nep->function = A;
676: }
677: if (B) {
678: PetscObjectReference((PetscObject)B);
679: MatDestroy(&nep->function_pre);
680: nep->function_pre = B;
681: }
682: nep->fui = NEP_USER_INTERFACE_CALLBACK;
683: nep->state = NEP_STATE_INITIAL;
684: return(0);
685: }
687: /*@C
688: NEPGetFunction - Returns the Function matrix and optionally the user
689: provided context for evaluating the Function.
691: Not Collective, but Mat object will be parallel if NEP object is
693: Input Parameter:
694: . nep - the nonlinear eigensolver context
696: Output Parameters:
697: + A - location to stash Function matrix (or NULL)
698: . B - location to stash preconditioner matrix (or NULL)
699: . fun - location to put Function function (or NULL)
700: - ctx - location to stash Function context (or NULL)
702: Level: advanced
704: .seealso: NEPSetFunction()
705: @*/
706: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,PetscErrorCode (**fun)(NEP,PetscScalar,Mat,Mat,void*),void **ctx)
707: {
710: NEPCheckCallback(nep,1);
711: if (A) *A = nep->function;
712: if (B) *B = nep->function_pre;
713: if (fun) *fun = nep->computefunction;
714: if (ctx) *ctx = nep->functionctx;
715: return(0);
716: }
718: /*@C
719: NEPSetJacobian - Sets the function to compute Jacobian T'(lambda) as well
720: as the location to store the matrix.
722: Logically Collective on NEP and Mat
724: Input Parameters:
725: + nep - the NEP context
726: . A - Jacobian matrix
727: . jac - Jacobian evaluation routine (if NULL then NEP retains any
728: previously set value)
729: - ctx - [optional] user-defined context for private data for the Jacobian
730: evaluation routine (may be NULL) (if NULL then NEP retains any
731: previously set value)
733: Calling Sequence of jac:
734: $ jac(NEP nep,PetscScalar lambda,Mat J,void *ctx)
736: + nep - the NEP context
737: . lambda - the scalar argument where T'(.) must be evaluated
738: . J - matrix that will contain T'(lambda)
739: - ctx - (optional) user-defined context, as set by NEPSetJacobian()
741: Level: beginner
743: .seealso: NEPSetFunction(), NEPGetJacobian()
744: @*/
745: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,PetscErrorCode (*jac)(NEP,PetscScalar,Mat,void*),void *ctx)
746: {
754: if (nep->state) { NEPReset(nep); }
755: else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }
757: if (jac) nep->computejacobian = jac;
758: if (ctx) nep->jacobianctx = ctx;
759: if (A) {
760: PetscObjectReference((PetscObject)A);
761: MatDestroy(&nep->jacobian);
762: nep->jacobian = A;
763: }
764: nep->fui = NEP_USER_INTERFACE_CALLBACK;
765: nep->state = NEP_STATE_INITIAL;
766: return(0);
767: }
769: /*@C
770: NEPGetJacobian - Returns the Jacobian matrix and optionally the user
771: provided routine and context for evaluating the Jacobian.
773: Not Collective, but Mat object will be parallel if NEP object is
775: Input Parameter:
776: . nep - the nonlinear eigensolver context
778: Output Parameters:
779: + A - location to stash Jacobian matrix (or NULL)
780: . jac - location to put Jacobian function (or NULL)
781: - ctx - location to stash Jacobian context (or NULL)
783: Level: advanced
785: .seealso: NEPSetJacobian()
786: @*/
787: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,PetscErrorCode (**jac)(NEP,PetscScalar,Mat,void*),void **ctx)
788: {
791: NEPCheckCallback(nep,1);
792: if (A) *A = nep->jacobian;
793: if (jac) *jac = nep->computejacobian;
794: if (ctx) *ctx = nep->jacobianctx;
795: return(0);
796: }
798: /*@
799: NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
800: in split form.
802: Collective on NEP, Mat and FN
804: Input Parameters:
805: + nep - the nonlinear eigensolver context
806: . n - number of terms in the split form
807: . A - array of matrices
808: . f - array of functions
809: - str - structure flag for matrices
811: Notes:
812: The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
813: for i=1,...,n. The derivative T'(lambda) can be obtained using the
814: derivatives of f_i.
816: The structure flag provides information about A_i's nonzero pattern
817: (see MatStructure enum). If all matrices have the same pattern, then
818: use SAME_NONZERO_PATTERN. If the patterns are different but contained
819: in the pattern of the first one, then use SUBSET_NONZERO_PATTERN.
820: Otherwise use DIFFERENT_NONZERO_PATTERN.
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()
828: @*/
829: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt n,Mat A[],FN f[],MatStructure str)
830: {
831: PetscInt i;
837: if (n <= 0) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %D",n);
843: for (i=0;i<n;i++) {
846: PetscObjectReference((PetscObject)A[i]);
847: PetscObjectReference((PetscObject)f[i]);
848: }
850: if (nep->state) { NEPReset(nep); }
851: else { NEPReset_Problem(nep); }
853: /* allocate space and copy matrices and functions */
854: PetscMalloc1(n,&nep->A);
855: PetscLogObjectMemory((PetscObject)nep,n*sizeof(Mat));
856: for (i=0;i<n;i++) nep->A[i] = A[i];
857: PetscMalloc1(n,&nep->f);
858: PetscLogObjectMemory((PetscObject)nep,n*sizeof(FN));
859: for (i=0;i<n;i++) nep->f[i] = f[i];
860: PetscCalloc1(n,&nep->nrma);
861: PetscLogObjectMemory((PetscObject)nep,n*sizeof(PetscReal));
862: nep->nt = n;
863: nep->mstr = str;
864: nep->fui = NEP_USER_INTERFACE_SPLIT;
865: nep->state = NEP_STATE_INITIAL;
866: return(0);
867: }
869: /*@
870: NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
871: the nonlinear operator in split form.
873: Not collective, though parallel Mats and FNs are returned if the NEP is parallel
875: Input Parameter:
876: + nep - the nonlinear eigensolver context
877: - k - the index of the requested term (starting in 0)
879: Output Parameters:
880: + A - the matrix of the requested term
881: - f - the function of the requested term
883: Level: intermediate
885: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
886: @*/
887: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
888: {
891: NEPCheckSplit(nep,1);
892: if (k<0 || k>=nep->nt) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",nep->nt-1);
893: if (A) *A = nep->A[k];
894: if (f) *f = nep->f[k];
895: return(0);
896: }
898: /*@
899: NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
900: the nonlinear operator, as well as the structure flag for matrices.
902: Not collective
904: Input Parameter:
905: . nep - the nonlinear eigensolver context
907: Output Parameters:
908: + n - the number of terms passed in NEPSetSplitOperator()
909: - str - the matrix structure flag passed in NEPSetSplitOperator()
911: Level: intermediate
913: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
914: @*/
915: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
916: {
919: NEPCheckSplit(nep,1);
920: if (n) *n = nep->nt;
921: if (str) *str = nep->mstr;
922: return(0);
923: }
925: /*@C
926: NEPSetDerivatives - Sets the function to compute the k-th derivative T^(k)(lambda)
927: for any value of k (including 0), as well as the location to store the matrix.
929: Logically Collective on NEP and Mat
931: Input Parameters:
932: + nep - the NEP context
933: . A - the matrix to store the computed derivative
934: . der - routing to evaluate the k-th derivative (if NULL then NEP retains any
935: previously set value)
936: - ctx - [optional] user-defined context for private data for the derivatives
937: evaluation routine (may be NULL) (if NULL then NEP retains any
938: previously set value)
940: Level: beginner
942: .seealso: NEPSetFunction(), NEPGetDerivatives()
943: @*/
944: PetscErrorCode NEPSetDerivatives(NEP nep,Mat A,PetscErrorCode (*der)(NEP,PetscScalar,PetscInt,Mat,void*),void *ctx)
945: {
953: if (nep->state) { NEPReset(nep); }
954: else { NEPReset_Problem(nep); }
956: if (der) nep->computederivatives = der;
957: if (ctx) nep->derivativesctx = ctx;
958: if (A) {
959: PetscObjectReference((PetscObject)A);
960: MatDestroy(&nep->derivatives);
961: nep->derivatives = A;
962: }
963: nep->fui = NEP_USER_INTERFACE_DERIVATIVES;
964: nep->state = NEP_STATE_INITIAL;
965: return(0);
966: }
968: /*@C
969: NEPGetDerivatives - Returns the derivatives matrix and optionally the user
970: provided routine and context for evaluating the derivatives.
972: Not Collective, but Mat object will be parallel if NEP object is
974: Input Parameter:
975: . nep - the nonlinear eigensolver context
977: Output Parameters:
978: + A - location to stash the derivatives matrix (or NULL)
979: . der - location to put derivatives function (or NULL)
980: - ctx - location to stash derivatives context (or NULL)
982: Level: advanced
984: .seealso: NEPSetDerivatives()
985: @*/
986: PetscErrorCode NEPGetDerivatives(NEP nep,Mat *A,PetscErrorCode (**der)(NEP,PetscScalar,PetscInt,Mat,void*),void **ctx)
987: {
990: NEPCheckDerivatives(nep,1);
991: if (A) *A = nep->derivatives;
992: if (der) *der = nep->computederivatives;
993: if (ctx) *ctx = nep->derivativesctx;
994: return(0);
995: }