Actual source code: interpol.c
slepc-3.21.1 2024-04-26
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: SLEPc nonlinear eigensolver: "interpol"
13: Method: Polynomial interpolation
15: Algorithm:
17: Uses a PEP object to solve the interpolated NEP. Currently supports
18: only Chebyshev interpolation on an interval.
20: References:
22: [1] C. Effenberger and D. Kresser, "Chebyshev interpolation for
23: nonlinear eigenvalue problems", BIT 52:933-951, 2012.
24: */
26: #include <slepc/private/nepimpl.h>
28: typedef struct {
29: PEP pep;
30: PetscReal tol; /* tolerance for norm of polynomial coefficients */
31: PetscInt maxdeg; /* maximum degree of interpolation polynomial */
32: PetscInt deg; /* actual degree of interpolation polynomial */
33: } NEP_INTERPOL;
35: static PetscErrorCode NEPSetUp_Interpol(NEP nep)
36: {
37: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
38: ST st;
39: RG rg;
40: PetscReal a,b,c,d,s,tol;
41: PetscScalar zero=0.0;
42: PetscBool flg,istrivial,trackall;
43: PetscInt its,in;
45: PetscFunctionBegin;
46: PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
47: PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
48: if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
49: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
50: PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
51: NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);
53: /* transfer PEP options */
54: if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
55: PetscCall(PEPSetBasis(ctx->pep,PEP_BASIS_CHEBYSHEV1));
56: PetscCall(PEPSetWhichEigenpairs(ctx->pep,PEP_TARGET_MAGNITUDE));
57: PetscCall(PetscObjectTypeCompare((PetscObject)ctx->pep,PEPJD,&flg));
58: if (!flg) {
59: PetscCall(PEPGetST(ctx->pep,&st));
60: PetscCall(STSetType(st,STSINVERT));
61: }
62: PetscCall(PEPSetDimensions(ctx->pep,nep->nev,nep->ncv,nep->mpd));
63: PetscCall(PEPGetTolerances(ctx->pep,&tol,&its));
64: if (tol==(PetscReal)PETSC_DEFAULT) tol = SlepcDefaultTol(nep->tol);
65: if (ctx->tol==(PetscReal)PETSC_DEFAULT) ctx->tol = tol;
66: if (its==PETSC_DEFAULT) its = nep->max_it;
67: PetscCall(PEPSetTolerances(ctx->pep,tol,its));
68: PetscCall(NEPGetTrackAll(nep,&trackall));
69: PetscCall(PEPSetTrackAll(ctx->pep,trackall));
71: /* transfer region options */
72: PetscCall(RGIsTrivial(nep->rg,&istrivial));
73: PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPINTERPOL requires a nontrivial region");
74: PetscCall(PetscObjectTypeCompare((PetscObject)nep->rg,RGINTERVAL,&flg));
75: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for interval regions");
76: PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,&c,&d));
77: PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for bounded intervals");
78: PetscCall(PEPGetRG(ctx->pep,&rg));
79: PetscCall(RGSetType(rg,RGINTERVAL));
80: PetscCheck(a!=b,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for intervals on the real axis");
81: s = 2.0/(b-a);
82: c = c*s;
83: d = d*s;
84: PetscCall(RGIntervalSetEndpoints(rg,-1.0,1.0,c,d));
85: PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
86: PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
87: PetscCall(PEPSetTarget(ctx->pep,(nep->target-(a+b)/2)*s));
89: PetscCall(NEPAllocateSolution(nep,0));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*
94: Input:
95: d, number of nodes to compute
96: a,b, interval extremes
97: Output:
98: *x, array containing the d Chebyshev nodes of the interval [a,b]
99: *dct2, coefficients to compute a discrete cosine transformation (DCT-II)
100: */
101: static PetscErrorCode ChebyshevNodes(PetscInt d,PetscReal a,PetscReal b,PetscScalar *x,PetscReal *dct2)
102: {
103: PetscInt j,i;
104: PetscReal t;
106: PetscFunctionBegin;
107: for (j=0;j<d+1;j++) {
108: t = ((2*j+1)*PETSC_PI)/(2*(d+1));
109: x[j] = (a+b)/2.0+((b-a)/2.0)*PetscCosReal(t);
110: for (i=0;i<d+1;i++) dct2[j*(d+1)+i] = PetscCosReal(i*t);
111: }
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode NEPSolve_Interpol(NEP nep)
116: {
117: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
118: Mat *A,*P;
119: PetscScalar *x,*fx,t;
120: PetscReal *cs,a,b,s,aprox,aprox0=1.0,*matnorm;
121: PetscInt i,j,k,deg=ctx->maxdeg;
122: PetscBool hasmnorm=PETSC_FALSE;
123: Vec vr,vi=NULL;
124: ST st;
126: PetscFunctionBegin;
127: PetscCall(PetscMalloc4((deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx,nep->nt,&matnorm));
128: for (j=0;j<nep->nt;j++) {
129: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
130: if (!hasmnorm) break;
131: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
132: }
133: if (!hasmnorm) for (j=0;j<nep->nt;j++) matnorm[j] = 1.0;
134: PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
135: PetscCall(ChebyshevNodes(deg,a,b,x,cs));
136: for (j=0;j<nep->nt;j++) {
137: for (i=0;i<=deg;i++) PetscCall(FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]));
138: }
139: /* Polynomial coefficients */
140: PetscCall(PetscMalloc1(deg+1,&A));
141: if (nep->P) PetscCall(PetscMalloc1(deg+1,&P));
142: ctx->deg = deg;
143: for (k=0;k<=deg;k++) {
144: PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]));
145: if (nep->P) PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P[k]));
146: t = 0.0;
147: for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
148: t *= 2.0/(deg+1);
149: if (k==0) t /= 2.0;
150: aprox = matnorm[0]*PetscAbsScalar(t);
151: PetscCall(MatScale(A[k],t));
152: if (nep->P) PetscCall(MatScale(P[k],t));
153: for (j=1;j<nep->nt;j++) {
154: t = 0.0;
155: for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
156: t *= 2.0/(deg+1);
157: if (k==0) t /= 2.0;
158: aprox += matnorm[j]*PetscAbsScalar(t);
159: PetscCall(MatAXPY(A[k],t,nep->A[j],nep->mstr));
160: if (nep->P) PetscCall(MatAXPY(P[k],t,nep->P[j],nep->mstrp));
161: }
162: if (k==0) aprox0 = aprox;
163: if (k>1 && aprox/aprox0<ctx->tol) { ctx->deg = k; deg = k; break; }
164: }
165: PetscCall(PEPSetOperators(ctx->pep,deg+1,A));
166: PetscCall(MatDestroyMatrices(deg+1,&A));
167: if (nep->P) {
168: PetscCall(PEPGetST(ctx->pep,&st));
169: PetscCall(STSetSplitPreconditioner(st,deg+1,P,nep->mstrp));
170: PetscCall(MatDestroyMatrices(deg+1,&P));
171: }
172: PetscCall(PetscFree4(cs,x,fx,matnorm));
174: /* Solve polynomial eigenproblem */
175: PetscCall(PEPSolve(ctx->pep));
176: PetscCall(PEPGetConverged(ctx->pep,&nep->nconv));
177: PetscCall(PEPGetIterationNumber(ctx->pep,&nep->its));
178: PetscCall(PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason));
179: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
180: PetscCall(BVCreateVec(nep->V,&vr));
181: #if !defined(PETSC_USE_COMPLEX)
182: PetscCall(VecDuplicate(vr,&vi));
183: #endif
184: s = 2.0/(b-a);
185: for (i=0;i<nep->nconv;i++) {
186: PetscCall(PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],vr,vi));
187: nep->eigr[i] /= s;
188: nep->eigr[i] += (a+b)/2.0;
189: nep->eigi[i] /= s;
190: PetscCall(BVInsertVec(nep->V,i,vr));
191: #if !defined(PETSC_USE_COMPLEX)
192: if (nep->eigi[i]!=0.0) PetscCall(BVInsertVec(nep->V,++i,vi));
193: #endif
194: }
195: PetscCall(VecDestroy(&vr));
196: PetscCall(VecDestroy(&vi));
198: nep->state = NEP_STATE_EIGENVECTORS;
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode PEPMonitor_Interpol(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
203: {
204: PetscInt i,n;
205: NEP nep = (NEP)ctx;
206: PetscReal a,b,s;
207: ST st;
209: PetscFunctionBegin;
210: n = PetscMin(nest,nep->ncv);
211: for (i=0;i<n;i++) {
212: nep->eigr[i] = eigr[i];
213: nep->eigi[i] = eigi[i];
214: nep->errest[i] = errest[i];
215: }
216: PetscCall(PEPGetST(pep,&st));
217: PetscCall(STBackTransform(st,n,nep->eigr,nep->eigi));
218: PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
219: s = 2.0/(b-a);
220: for (i=0;i<n;i++) {
221: nep->eigr[i] /= s;
222: nep->eigr[i] += (a+b)/2.0;
223: nep->eigi[i] /= s;
224: }
225: PetscCall(NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode NEPSetFromOptions_Interpol(NEP nep,PetscOptionItems *PetscOptionsObject)
230: {
231: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
232: PetscInt i;
233: PetscBool flg1,flg2;
234: PetscReal r;
236: PetscFunctionBegin;
237: PetscOptionsHeadBegin(PetscOptionsObject,"NEP Interpol Options");
239: PetscCall(NEPInterpolGetInterpolation(nep,&r,&i));
240: if (!i) i = PETSC_DEFAULT;
241: PetscCall(PetscOptionsInt("-nep_interpol_interpolation_degree","Maximum degree of polynomial interpolation","NEPInterpolSetInterpolation",i,&i,&flg1));
242: PetscCall(PetscOptionsReal("-nep_interpol_interpolation_tol","Tolerance for interpolation coefficients","NEPInterpolSetInterpolation",r,&r,&flg2));
243: if (flg1 || flg2) PetscCall(NEPInterpolSetInterpolation(nep,r,i));
245: PetscOptionsHeadEnd();
247: if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
248: PetscCall(PEPSetFromOptions(ctx->pep));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode NEPInterpolSetInterpolation_Interpol(NEP nep,PetscReal tol,PetscInt degree)
253: {
254: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
256: PetscFunctionBegin;
257: if (tol == (PetscReal)PETSC_DEFAULT) {
258: ctx->tol = PETSC_DEFAULT;
259: nep->state = NEP_STATE_INITIAL;
260: } else {
261: PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
262: ctx->tol = tol;
263: }
264: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
265: ctx->maxdeg = 0;
266: if (nep->state) PetscCall(NEPReset(nep));
267: nep->state = NEP_STATE_INITIAL;
268: } else {
269: PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
270: if (ctx->maxdeg != degree) {
271: ctx->maxdeg = degree;
272: if (nep->state) PetscCall(NEPReset(nep));
273: nep->state = NEP_STATE_INITIAL;
274: }
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: NEPInterpolSetInterpolation - Sets the tolerance and maximum degree when building
281: the interpolation polynomial.
283: Collective
285: Input Parameters:
286: + nep - nonlinear eigenvalue solver
287: . tol - tolerance to stop computing polynomial coefficients
288: - deg - maximum degree of interpolation
290: Options Database Key:
291: + -nep_interpol_interpolation_tol <tol> - Sets the tolerance to stop computing polynomial coefficients
292: - -nep_interpol_interpolation_degree <degree> - Sets the maximum degree of interpolation
294: Notes:
295: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
297: Level: advanced
299: .seealso: NEPInterpolGetInterpolation()
300: @*/
301: PetscErrorCode NEPInterpolSetInterpolation(NEP nep,PetscReal tol,PetscInt deg)
302: {
303: PetscFunctionBegin;
307: PetscTryMethod(nep,"NEPInterpolSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,deg));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode NEPInterpolGetInterpolation_Interpol(NEP nep,PetscReal *tol,PetscInt *deg)
312: {
313: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
315: PetscFunctionBegin;
316: if (tol) *tol = ctx->tol;
317: if (deg) *deg = ctx->maxdeg;
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@
322: NEPInterpolGetInterpolation - Gets the tolerance and maximum degree when building
323: the interpolation polynomial.
325: Not Collective
327: Input Parameter:
328: . nep - nonlinear eigenvalue solver
330: Output Parameters:
331: + tol - tolerance to stop computing polynomial coefficients
332: - deg - maximum degree of interpolation
334: Level: advanced
336: .seealso: NEPInterpolSetInterpolation()
337: @*/
338: PetscErrorCode NEPInterpolGetInterpolation(NEP nep,PetscReal *tol,PetscInt *deg)
339: {
340: PetscFunctionBegin;
342: PetscUseMethod(nep,"NEPInterpolGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,deg));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode NEPInterpolSetPEP_Interpol(NEP nep,PEP pep)
347: {
348: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
350: PetscFunctionBegin;
351: PetscCall(PetscObjectReference((PetscObject)pep));
352: PetscCall(PEPDestroy(&ctx->pep));
353: ctx->pep = pep;
354: nep->state = NEP_STATE_INITIAL;
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: NEPInterpolSetPEP - Associate a polynomial eigensolver object (PEP) to the
360: nonlinear eigenvalue solver.
362: Collective
364: Input Parameters:
365: + nep - nonlinear eigenvalue solver
366: - pep - the polynomial eigensolver object
368: Level: advanced
370: .seealso: NEPInterpolGetPEP()
371: @*/
372: PetscErrorCode NEPInterpolSetPEP(NEP nep,PEP pep)
373: {
374: PetscFunctionBegin;
377: PetscCheckSameComm(nep,1,pep,2);
378: PetscTryMethod(nep,"NEPInterpolSetPEP_C",(NEP,PEP),(nep,pep));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode NEPInterpolGetPEP_Interpol(NEP nep,PEP *pep)
383: {
384: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
386: PetscFunctionBegin;
387: if (!ctx->pep) {
388: PetscCall(PEPCreate(PetscObjectComm((PetscObject)nep),&ctx->pep));
389: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->pep,(PetscObject)nep,1));
390: PetscCall(PEPSetOptionsPrefix(ctx->pep,((PetscObject)nep)->prefix));
391: PetscCall(PEPAppendOptionsPrefix(ctx->pep,"nep_interpol_"));
392: PetscCall(PetscObjectSetOptions((PetscObject)ctx->pep,((PetscObject)nep)->options));
393: PetscCall(PEPMonitorSet(ctx->pep,PEPMonitor_Interpol,nep,NULL));
394: }
395: *pep = ctx->pep;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /*@
400: NEPInterpolGetPEP - Retrieve the polynomial eigensolver object (PEP)
401: associated with the nonlinear eigenvalue solver.
403: Collective
405: Input Parameter:
406: . nep - nonlinear eigenvalue solver
408: Output Parameter:
409: . pep - the polynomial eigensolver object
411: Level: advanced
413: .seealso: NEPInterpolSetPEP()
414: @*/
415: PetscErrorCode NEPInterpolGetPEP(NEP nep,PEP *pep)
416: {
417: PetscFunctionBegin;
419: PetscAssertPointer(pep,2);
420: PetscUseMethod(nep,"NEPInterpolGetPEP_C",(NEP,PEP*),(nep,pep));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: static PetscErrorCode NEPView_Interpol(NEP nep,PetscViewer viewer)
425: {
426: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
427: PetscBool isascii;
429: PetscFunctionBegin;
430: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
431: if (isascii) {
432: if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
433: PetscCall(PetscViewerASCIIPrintf(viewer," polynomial degree %" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->deg,ctx->maxdeg));
434: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance for norm of polynomial coefficients %g\n",(double)ctx->tol));
435: PetscCall(PetscViewerASCIIPushTab(viewer));
436: PetscCall(PEPView(ctx->pep,viewer));
437: PetscCall(PetscViewerASCIIPopTab(viewer));
438: }
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: static PetscErrorCode NEPReset_Interpol(NEP nep)
443: {
444: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
446: PetscFunctionBegin;
447: PetscCall(PEPReset(ctx->pep));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode NEPDestroy_Interpol(NEP nep)
452: {
453: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
455: PetscFunctionBegin;
456: PetscCall(PEPDestroy(&ctx->pep));
457: PetscCall(PetscFree(nep->data));
458: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NULL));
459: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NULL));
460: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NULL));
461: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NULL));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: SLEPC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
466: {
467: NEP_INTERPOL *ctx;
469: PetscFunctionBegin;
470: PetscCall(PetscNew(&ctx));
471: nep->data = (void*)ctx;
472: ctx->maxdeg = 5;
473: ctx->tol = PETSC_DEFAULT;
475: nep->ops->solve = NEPSolve_Interpol;
476: nep->ops->setup = NEPSetUp_Interpol;
477: nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
478: nep->ops->reset = NEPReset_Interpol;
479: nep->ops->destroy = NEPDestroy_Interpol;
480: nep->ops->view = NEPView_Interpol;
482: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NEPInterpolSetInterpolation_Interpol));
483: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NEPInterpolGetInterpolation_Interpol));
484: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol));
485: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }