Actual source code: stoar.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: SLEPc polynomial eigensolver: "stoar"
13: Method: S-TOAR
15: Algorithm:
17: Symmetric Two-Level Orthogonal Arnoldi.
19: References:
21: [1] C. Campos and J.E. Roman, "Restarted Q-Arnoldi-type methods
22: exploiting symmetry in quadratic eigenvalue problems", BIT
23: Numer. Math. 56(4):1213-1236, 2016.
24: */
26: #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
27: #include "../src/pep/impls/krylov/pepkrylov.h"
28: #include <slepcblaslapack.h>
30: static PetscBool cited = PETSC_FALSE;
31: static const char citation[] =
32: "@Article{slepc-stoar,\n"
33: " author = \"C. Campos and J. E. Roman\",\n"
34: " title = \"Restarted {Q-Arnoldi-type} methods exploiting symmetry in quadratic eigenvalue problems\",\n"
35: " journal = \"{BIT} Numer. Math.\",\n"
36: " volume = \"56\",\n"
37: " number = \"4\",\n"
38: " pages = \"1213--1236\",\n"
39: " year = \"2016,\"\n"
40: " doi = \"https://doi.org/10.1007/s10543-016-0601-5\"\n"
41: "}\n";
43: typedef struct {
44: PetscReal scal[2];
45: Mat A[2];
46: Vec t;
47: } ShellMatCtx;
49: PetscErrorCode MatMult_Func(Mat A,Vec x,Vec y)
50: {
52: ShellMatCtx *ctx;
55: MatShellGetContext(A,&ctx);
56: MatMult(ctx->A[0],x,y);
57: VecScale(y,ctx->scal[0]);
58: if (ctx->scal[1]) {
59: MatMult(ctx->A[1],x,ctx->t);
60: VecAXPY(y,ctx->scal[1],ctx->t);
61: }
62: return(0);
63: }
65: PetscErrorCode MatDestroy_Func(Mat A)
66: {
67: ShellMatCtx *ctx;
71: MatShellGetContext(A,(void**)&ctx);
72: VecDestroy(&ctx->t);
73: PetscFree(ctx);
74: return(0);
75: }
77: PetscErrorCode PEPSTOARSetUpInnerMatrix(PEP pep,Mat *B)
78: {
79: Mat pB[4],Bs[3],D[3];
80: PetscInt i,j,n,m;
81: ShellMatCtx *ctxMat[3];
82: PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
86: for (i=0;i<3;i++) {
87: STGetMatrixTransformed(pep->st,i,&D[i]); /* D[2] = M */
88: }
89: MatGetLocalSize(D[2],&m,&n);
90:
91: for (j=0;j<3;j++) {
92: PetscNew(ctxMat+j);
93: MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&Bs[j]);
94: MatShellSetOperation(Bs[j],MATOP_MULT,(void(*)(void))MatMult_Func);
95: MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_Func);
96: }
97: for (i=0;i<4;i++) pB[i] = NULL;
98: if (ctx->alpha) {
99: ctxMat[0]->A[0] = D[0]; ctxMat[0]->scal[0] = ctx->alpha; ctxMat[0]->scal[1] = 0.0;
100: ctxMat[2]->A[0] = D[2]; ctxMat[2]->scal[0] = -ctx->alpha*pep->sfactor*pep->sfactor; ctxMat[2]->scal[1] = 0.0;
101: pB[0] = Bs[0]; pB[3] = Bs[2];
102: }
103: if (ctx->beta) {
104: i = (ctx->alpha)?1:0;
105: ctxMat[0]->scal[1] = 0.0;
106: ctxMat[0]->A[i] = D[1]; ctxMat[0]->scal[i] = -ctx->beta*pep->sfactor;
107: ctxMat[1]->A[0] = D[2]; ctxMat[1]->scal[0] = -ctx->beta*pep->sfactor*pep->sfactor; ctxMat[1]->scal[1] = 0.0;
108: pB[0] = Bs[0]; pB[1] = pB[2] = Bs[1];
109: }
110: BVCreateVec(pep->V,&ctxMat[0]->t);
111: MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pB,B);
112: for (j=0;j<3;j++) { MatDestroy(&Bs[j]); }
113: return(0);
114: }
116: PetscErrorCode PEPSetUp_STOAR(PEP pep)
117: {
118: PetscErrorCode ierr;
119: PetscBool shift,sinv,flg;
120: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
121: PetscInt ld,i;
122: PetscReal eta;
123: BVOrthogType otype;
124: BVOrthogBlockType obtype;
127: if (pep->problem_type!=PEP_HERMITIAN && pep->problem_type!=PEP_HYPERBOLIC) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Requested method is only available for Hermitian (or hyperbolic) problems");
128: if (pep->nmat!=3) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver only available for quadratic problems");
129: if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver not implemented for non-monomial bases");
130: /* spectrum slicing requires special treatment of default values */
131: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
132: if (pep->which==PEP_ALL) {
133: if (pep->stopping!=PEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Spectrum slicing does not support user-defined stopping test");
134: pep->ops->solve = PEPSolve_STOAR_QSlice;
135: pep->ops->extractvectors = NULL;
136: pep->ops->setdefaultst = NULL;
137: PEPSetUp_STOAR_QSlice(pep);
138: } else {
139: PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
140: if (!ctx->lock && pep->mpd<pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
141: if (!pep->max_it) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
142: pep->ops->solve = PEPSolve_STOAR;
143: ld = pep->ncv+2;
144: DSSetType(pep->ds,DSGHIEP);
145: DSSetCompact(pep->ds,PETSC_TRUE);
146: DSAllocate(pep->ds,ld);
147: PEPBasisCoefficients(pep,pep->pbc);
148: STGetTransform(pep->st,&flg);
149: if (!flg) {
150: PetscFree(pep->solvematcoeffs);
151: PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
152: PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
153: if (sinv) {
154: PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
155: } else {
156: for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
157: pep->solvematcoeffs[pep->nmat-1] = 1.0;
158: }
159: }
160: }
162: pep->lineariz = PETSC_TRUE;
163: PetscObjectTypeCompare((PetscObject)pep->st,STSHIFT,&shift);
164: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
165: if (!shift && !sinv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only STSHIFT and STSINVERT spectral transformations can be used");
166: if (!pep->which) {
167: if (sinv) pep->which = PEP_TARGET_MAGNITUDE;
168: else pep->which = PEP_LARGEST_MAGNITUDE;
169: }
171: PEPAllocateSolution(pep,2);
172: PEPSetWorkVecs(pep,4);
173: BVDestroy(&ctx->V);
174: BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
175: BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype);
176: BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
177: return(0);
178: }
180: /*
181: Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
182: */
183: static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
184: {
186: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
187: PetscInt i,j,m=*M,l,lock;
188: PetscInt lds,d,ld,offq,nqt,ldds;
189: Vec v=t_[0],t=t_[1],q=t_[2];
190: PetscReal norm,sym=0.0,fro=0.0,*f;
191: PetscScalar *y,*S,*x,sigma;
192: PetscBLASInt j_,one=1;
193: PetscBool lindep,flg,sinvert=PETSC_FALSE;
194: Mat MS;
197: PetscMalloc1(*M,&y);
198: BVGetSizes(pep->V,NULL,NULL,&ld);
199: BVTensorGetDegree(ctx->V,&d);
200: BVGetActiveColumns(pep->V,&lock,&nqt);
201: lds = d*ld;
202: offq = ld;
203: DSGetLeadingDimension(pep->ds,&ldds);
204: *breakdown = PETSC_FALSE; /* ----- */
205: DSGetDimensions(pep->ds,NULL,NULL,&l,NULL,NULL);
206: BVSetActiveColumns(ctx->V,0,m);
207: BVSetActiveColumns(pep->V,0,nqt);
208: STGetTransform(pep->st,&flg);
209: if (!flg) {
210: /* spectral transformation handled by the solver */
211: PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,"");
212: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"ST type not supported for TOAR without transforming matrices");
213: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert);
214: STGetShift(pep->st,&sigma);
215: }
216: for (j=k;j<m;j++) {
217: /* apply operator */
218: BVTensorGetFactors(ctx->V,NULL,&MS);
219: MatDenseGetArray(MS,&S);
220: BVGetColumn(pep->V,nqt,&t);
221: BVMultVec(pep->V,1.0,0.0,v,S+j*lds);
222: if (!sinvert) {
223: STMatMult(pep->st,0,v,q);
224: BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
225: STMatMult(pep->st,1,v,t);
226: VecAXPY(q,pep->sfactor,t);
227: if (ctx->beta && ctx->alpha) {
228: STMatMult(pep->st,2,v,t);
229: VecAXPY(q,-pep->sfactor*pep->sfactor*ctx->beta/ctx->alpha,t);
230: }
231: STMatSolve(pep->st,q,t);
232: VecScale(t,-1.0/(pep->sfactor*pep->sfactor));
233: } else {
234: STMatMult(pep->st,1,v,q);
235: STMatMult(pep->st,2,v,t);
236: VecAXPY(q,sigma*pep->sfactor,t);
237: VecScale(q,pep->sfactor);
238: BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
239: STMatMult(pep->st,2,v,t);
240: VecAXPY(q,pep->sfactor*pep->sfactor,t);
241: STMatSolve(pep->st,q,t);
242: VecScale(t,-1.0);
243: }
244: BVRestoreColumn(pep->V,nqt,&t);
246: /* orthogonalize */
247: if (!sinvert) x = S+offq+(j+1)*lds;
248: else x = S+(j+1)*lds;
249: BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep);
251: if (!lindep) {
252: if (!sinvert) *(S+offq+(j+1)*lds+nqt) = norm;
253: else *(S+(j+1)*lds+nqt) = norm;
254: BVScaleColumn(pep->V,nqt,1.0/norm);
255: nqt++;
256: }
257: if (!sinvert) {
258: for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
259: if (ctx->beta && ctx->alpha) {
260: for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+offq+i) -= *(S+(j+1)*lds+i)*ctx->beta/ctx->alpha;
261: }
262: } else for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
263: BVSetActiveColumns(pep->V,0,nqt);
264: MatDenseRestoreArray(MS,&S);
265: BVTensorRestoreFactors(ctx->V,NULL,&MS);
267: /* level-2 orthogonalization */
268: BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep);
269: a[j] = PetscRealPart(y[j]);
270: omega[j+1] = (norm > 0)?1.0:-1.0;
271: BVScaleColumn(ctx->V,j+1,1.0/norm);
272: b[j] = PetscAbsReal(norm);
274: /* check symmetry */
275: DSGetArrayReal(pep->ds,DS_MAT_T,&f);
276: if (j==k) {
277: for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
278: for (i=0;i<l;i++) y[i] = 0.0;
279: }
280: DSRestoreArrayReal(pep->ds,DS_MAT_T,&f);
281: if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
282: PetscBLASIntCast(j,&j_);
283: sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
284: fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
285: if (j>0) fro = SlepcAbs(fro,b[j-1]);
286: if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
287: *symmlost = PETSC_TRUE;
288: *M=j;
289: break;
290: }
291: }
292: BVSetActiveColumns(pep->V,lock,nqt);
293: BVSetActiveColumns(ctx->V,0,*M);
294: PetscFree(y);
295: return(0);
296: }
298: #if 0
299: static PetscErrorCode PEPSTOARpreKConvergence(PEP pep,PetscInt nv,PetscReal *norm,Vec *w)
300: {
302: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
303: PetscBLASInt n_,one=1;
304: PetscInt lds=2*ctx->ld;
305: PetscReal t1,t2;
306: PetscScalar *S=ctx->S;
309: PetscBLASIntCast(nv+2,&n_);
310: t1 = BLASnrm2_(&n_,S+nv*2*ctx->ld,&one);
311: t2 = BLASnrm2_(&n_,S+(nv*2+1)*ctx->ld,&one);
312: *norm = SlepcAbs(t1,t2);
313: BVSetActiveColumns(pep->V,0,nv+2);
314: BVMultVec(pep->V,1.0,0.0,w[1],S+nv*lds);
315: STMatMult(pep->st,0,w[1],w[2]);
316: VecNorm(w[2],NORM_2,&t1);
317: BVMultVec(pep->V,1.0,0.0,w[1],S+ctx->ld+nv*lds);
318: STMatMult(pep->st,2,w[1],w[2]);
319: VecNorm(w[2],NORM_2,&t2);
320: t2 *= pep->sfactor*pep->sfactor;
321: *norm = PetscMax(*norm,SlepcAbs(t1,t2));
322: return(0);
323: }
324: #endif
326: PetscErrorCode PEPSolve_STOAR(PEP pep)
327: {
329: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
330: PetscInt j,k,l,nv=0,ld,ldds,t,nq=0;
331: PetscInt nconv=0,deg=pep->nmat-1;
332: PetscScalar *Q,*om,sigma;
333: PetscReal beta,norm=1.0,*omega,*a,*b,*r;
334: PetscBool breakdown,symmlost=PETSC_FALSE,sinv,falselock=PETSC_TRUE,flg;
335: Mat MQ,A;
336: Vec vomega;
339: PetscCitationsRegister(citation,&cited);
340: PEPSTOARSetUpInnerMatrix(pep,&A);
341: BVSetMatrix(ctx->V,A,PETSC_TRUE);
342: MatDestroy(&A);
343: if (ctx->lock) {
344: PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL);
345: }
346: BVGetSizes(pep->V,NULL,NULL,&ld);
347: STGetShift(pep->st,&sigma);
348: STGetTransform(pep->st,&flg);
349: if (pep->sfactor!=1.0) {
350: if (!flg) {
351: pep->target /= pep->sfactor;
352: RGPushScale(pep->rg,1.0/pep->sfactor);
353: STScaleShift(pep->st,1.0/pep->sfactor);
354: sigma /= pep->sfactor;
355: } else {
356: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
357: pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
358: RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
359: STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
360: }
361: }
362: if (flg) sigma = 0.0;
364: /* Get the starting Arnoldi vector */
365: BVTensorBuildFirstColumn(ctx->V,pep->nini);
366: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
367: VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
368: BVSetActiveColumns(ctx->V,0,1);
369: BVGetSignature(ctx->V,vomega);
370: VecGetArray(vomega,&om);
371: omega[0] = PetscRealPart(om[0]);
372: VecRestoreArray(vomega,&om);
373: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
374: VecDestroy(&vomega);
376: /* Restart loop */
377: l = 0;
378: DSGetLeadingDimension(pep->ds,&ldds);
379: while (pep->reason == PEP_CONVERGED_ITERATING) {
380: pep->its++;
381: DSGetArrayReal(pep->ds,DS_MAT_T,&a);
382: b = a+ldds;
383: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
385: /* Compute an nv-step Lanczos factorization */
386: nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
387: PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work);
388: beta = b[nv-1];
389: if (symmlost && nv==pep->nconv+l) {
390: pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
391: pep->nconv = nconv;
392: if (falselock || !ctx->lock) {
393: BVSetActiveColumns(ctx->V,0,pep->nconv);
394: BVTensorCompress(ctx->V,0);
395: }
396: break;
397: }
398: DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
399: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
400: DSSetDimensions(pep->ds,nv,0,pep->nconv,pep->nconv+l);
401: if (l==0) {
402: DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
403: } else {
404: DSSetState(pep->ds,DS_STATE_RAW);
405: }
407: /* Solve projected problem */
408: DSSolve(pep->ds,pep->eigr,pep->eigi);
409: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
410: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
412: /* Check convergence */
413: /* PEPSTOARpreKConvergence(pep,nv,&norm,pep->work);*/
414: norm = 1.0;
415: DSGetDimensions(pep->ds,NULL,NULL,NULL,NULL,&t);
416: PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k);
417: (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
419: /* Update l */
420: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
421: else {
422: l = PetscMax(1,(PetscInt)((nv-k)/2));
423: l = PetscMin(l,t);
424: if (!breakdown) {
425: DSGetArrayReal(pep->ds,DS_MAT_T,&a);
426: if (*(a+ldds+k+l-1)!=0) {
427: if (k+l<nv-1) l = l+1;
428: else l = l-1;
429: }
430: /* Prepare the Rayleigh quotient for restart */
431: DSGetArray(pep->ds,DS_MAT_Q,&Q);
432: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
433: r = a + 2*ldds;
434: for (j=k;j<k+l;j++) {
435: r[j] = PetscRealPart(Q[nv-1+j*ldds]*beta);
436: }
437: b = a+ldds;
438: b[k+l-1] = r[k+l-1];
439: omega[k+l] = omega[nv];
440: DSRestoreArray(pep->ds,DS_MAT_Q,&Q);
441: DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
442: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
443: }
444: }
445: nconv = k;
446: if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
448: /* Update S */
449: DSGetMat(pep->ds,DS_MAT_Q,&MQ);
450: BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
451: MatDestroy(&MQ);
453: /* Copy last column of S */
454: BVCopyColumn(ctx->V,nv,k+l);
455: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
456: VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
457: VecGetArray(vomega,&om);
458: for (j=0;j<k+l;j++) om[j] = omega[j];
459: VecRestoreArray(vomega,&om);
460: BVSetActiveColumns(ctx->V,0,k+l);
461: BVSetSignature(ctx->V,vomega);
462: VecDestroy(&vomega);
463: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
465: if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
466: /* stop if breakdown */
467: PetscInfo2(pep,"Breakdown TOAR method (it=%D norm=%g)\n",pep->its,(double)beta);
468: pep->reason = PEP_DIVERGED_BREAKDOWN;
469: }
470: if (pep->reason != PEP_CONVERGED_ITERATING) l--;
471: BVGetActiveColumns(pep->V,NULL,&nq);
472: if (k+l+deg<=nq) {
473: BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
474: if (!falselock && ctx->lock) {
475: BVTensorCompress(ctx->V,k-pep->nconv);
476: } else {
477: BVTensorCompress(ctx->V,0);
478: }
479: }
480: pep->nconv = k;
481: PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
482: }
484: if (pep->nconv>0) {
485: BVSetActiveColumns(ctx->V,0,pep->nconv);
486: BVGetActiveColumns(pep->V,NULL,&nq);
487: BVSetActiveColumns(pep->V,0,nq);
488: if (nq>pep->nconv) {
489: BVTensorCompress(ctx->V,pep->nconv);
490: BVSetActiveColumns(pep->V,0,pep->nconv);
491: }
492: }
493: STGetTransform(pep->st,&flg);
494: if (!flg && pep->ops->backtransform) {
495: (*pep->ops->backtransform)(pep);
496: }
497: if (pep->sfactor!=1.0) {
498: for (j=0;j<pep->nconv;j++) {
499: pep->eigr[j] *= pep->sfactor;
500: pep->eigi[j] *= pep->sfactor;
501: }
502: }
503: /* restore original values */
504: if (!flg) {
505: pep->target *= pep->sfactor;
506: STScaleShift(pep->st,pep->sfactor);
507: } else {
508: STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
509: pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
510: }
511: if (pep->sfactor!=1.0) { RGPopScale(pep->rg); }
513: /* truncate Schur decomposition and change the state to raw so that
514: DSVectors() computes eigenvectors from scratch */
515: DSSetDimensions(pep->ds,pep->nconv,0,0,0);
516: DSSetState(pep->ds,DS_STATE_RAW);
517: return(0);
518: }
520: PetscErrorCode PEPSetFromOptions_STOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
521: {
523: PetscBool flg,lock,b,f1,f2,f3;
524: PetscInt i,j,k;
525: PetscReal array[2]={0,0};
526: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
529: PetscOptionsHead(PetscOptionsObject,"PEP STOAR Options");
531: PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg);
532: if (flg) { PEPSTOARSetLocking(pep,lock); }
534: b = ctx->detect;
535: PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg);
536: if (flg) { PEPSTOARSetDetectZeros(pep,b); }
538: i = 1;
539: j = k = PETSC_DECIDE;
540: PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1);
541: PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2);
542: PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3);
543: if (f1 || f2 || f3) { PEPSTOARSetDimensions(pep,i,j,k); }
545: k = 2;
546: PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg);
547: if (flg) {
548: PEPSTOARSetLinearization(pep,array[0],array[1]);
549: }
551: b = ctx->checket;
552: PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg);
553: if (flg) { PEPSTOARSetCheckEigenvalueType(pep,b); }
555: PetscOptionsTail();
556: return(0);
557: }
559: static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
560: {
561: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
564: ctx->lock = lock;
565: return(0);
566: }
568: /*@
569: PEPSTOARSetLocking - Choose between locking and non-locking variants of
570: the STOAR method.
572: Logically Collective on PEP
574: Input Parameters:
575: + pep - the eigenproblem solver context
576: - lock - true if the locking variant must be selected
578: Options Database Key:
579: . -pep_stoar_locking - Sets the locking flag
581: Notes:
582: The default is to lock converged eigenpairs when the method restarts.
583: This behaviour can be changed so that all directions are kept in the
584: working subspace even if already converged to working accuracy (the
585: non-locking variant).
587: Level: advanced
589: .seealso: PEPSTOARGetLocking()
590: @*/
591: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
592: {
598: PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
599: return(0);
600: }
602: static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
603: {
604: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
607: *lock = ctx->lock;
608: return(0);
609: }
611: /*@
612: PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.
614: Not Collective
616: Input Parameter:
617: . pep - the eigenproblem solver context
619: Output Parameter:
620: . lock - the locking flag
622: Level: advanced
624: .seealso: PEPSTOARSetLocking()
625: @*/
626: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
627: {
633: PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
634: return(0);
635: }
637: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
638: {
640: PetscInt i,numsh;
641: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
642: PEP_SR sr = ctx->sr;
645: if (!pep->state) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
646: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
647: switch (pep->state) {
648: case PEP_STATE_INITIAL:
649: break;
650: case PEP_STATE_SETUP:
651: if (n) *n = 2;
652: if (shifts) {
653: PetscMalloc1(2,shifts);
654: (*shifts)[0] = pep->inta;
655: (*shifts)[1] = pep->intb;
656: }
657: if (inertias) {
658: PetscMalloc1(2,inertias);
659: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
660: (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
661: }
662: break;
663: case PEP_STATE_SOLVED:
664: case PEP_STATE_EIGENVECTORS:
665: numsh = ctx->nshifts;
666: if (n) *n = numsh;
667: if (shifts) {
668: PetscMalloc1(numsh,shifts);
669: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
670: }
671: if (inertias) {
672: PetscMalloc1(numsh,inertias);
673: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
674: }
675: break;
676: }
677: return(0);
678: }
680: /*@C
681: PEPSTOARGetInertias - Gets the values of the shifts and their
682: corresponding inertias in case of doing spectrum slicing for a
683: computational interval.
685: Not Collective
687: Input Parameter:
688: . pep - the eigenproblem solver context
690: Output Parameters:
691: + n - number of shifts, including the endpoints of the interval
692: . shifts - the values of the shifts used internally in the solver
693: - inertias - the values of the inertia in each shift
695: Notes:
696: If called after PEPSolve(), all shifts used internally by the solver are
697: returned (including both endpoints and any intermediate ones). If called
698: before PEPSolve() and after PEPSetUp() then only the information of the
699: endpoints of subintervals is available.
701: This function is only available for spectrum slicing runs.
703: The returned arrays should be freed by the user. Can pass NULL in any of
704: the two arrays if not required.
706: Fortran Notes:
707: The calling sequence from Fortran is
708: .vb
709: PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
710: integer n
711: double precision shifts(*)
712: integer inertias(*)
713: .ve
714: The arrays should be at least of length n. The value of n can be determined
715: by an initial call
716: .vb
717: PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
718: .ve
720: Level: advanced
722: .seealso: PEPSetInterval()
723: @*/
724: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
725: {
731: PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
732: return(0);
733: }
735: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
736: {
737: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
740: ctx->detect = detect;
741: pep->state = PEP_STATE_INITIAL;
742: return(0);
743: }
745: /*@
746: PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
747: zeros during the factorizations throughout the spectrum slicing computation.
749: Logically Collective on PEP
751: Input Parameters:
752: + pep - the eigenproblem solver context
753: - detect - check for zeros
755: Options Database Key:
756: . -pep_stoar_detect_zeros - Check for zeros; this takes an optional
757: bool value (0/1/no/yes/true/false)
759: Notes:
760: A zero in the factorization indicates that a shift coincides with an eigenvalue.
762: This flag is turned off by default, and may be necessary in some cases.
763: This feature currently requires an external package for factorizations
764: with support for zero detection, e.g. MUMPS.
766: Level: advanced
768: .seealso: PEPSetInterval()
769: @*/
770: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
771: {
777: PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
778: return(0);
779: }
781: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
782: {
783: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
786: *detect = ctx->detect;
787: return(0);
788: }
790: /*@
791: PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
792: in spectrum slicing.
794: Not Collective
796: Input Parameter:
797: . pep - the eigenproblem solver context
799: Output Parameter:
800: . detect - whether zeros detection is enforced during factorizations
802: Level: advanced
804: .seealso: PEPSTOARSetDetectZeros()
805: @*/
806: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
807: {
813: PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
814: return(0);
815: }
817: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
818: {
819: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
822: if (beta==0.0 && alpha==0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
823: ctx->alpha = alpha;
824: ctx->beta = beta;
825: return(0);
826: }
828: /*@
829: PEPSTOARSetLinearization - Set the coefficients that define
830: the linearization of a quadratic eigenproblem.
832: Logically Collective on PEP
834: Input Parameters:
835: + pep - polynomial eigenvalue solver
836: . alpha - first parameter of the linearization
837: - beta - second parameter of the linearization
839: Options Database Key:
840: . -pep_stoar_linearization <alpha,beta> - Sets the coefficients
842: Notes:
843: Cannot pass zero for both alpha and beta. The default values are
844: alpha=1 and beta=0.
846: Level: advanced
848: .seealso: PEPSTOARGetLinearization()
849: @*/
850: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
851: {
858: PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
859: return(0);
860: }
862: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
863: {
864: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
867: if (alpha) *alpha = ctx->alpha;
868: if (beta) *beta = ctx->beta;
869: return(0);
870: }
872: /*@
873: PEPSTOARGetLinearization - Returns the coefficients that define
874: the linearization of a quadratic eigenproblem.
876: Not Collective
878: Input Parameter:
879: . pep - polynomial eigenvalue solver
881: Output Parameters:
882: + alpha - the first parameter of the linearization
883: - beta - the second parameter of the linearization
885: Level: advanced
887: .seealso: PEPSTOARSetLinearization()
888: @*/
889: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
890: {
895: PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
896: return(0);
897: }
899: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
900: {
901: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
904: if (nev<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
905: ctx->nev = nev;
906: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
907: ctx->ncv = 0;
908: } else {
909: if (ncv<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
910: ctx->ncv = ncv;
911: }
912: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
913: ctx->mpd = 0;
914: } else {
915: if (mpd<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
916: ctx->mpd = mpd;
917: }
918: pep->state = PEP_STATE_INITIAL;
919: return(0);
920: }
922: /*@
923: PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
924: step in case of doing spectrum slicing for a computational interval.
925: The meaning of the parameters is the same as in PEPSetDimensions().
927: Logically Collective on PEP
929: Input Parameters:
930: + pep - the eigenproblem solver context
931: . nev - number of eigenvalues to compute
932: . ncv - the maximum dimension of the subspace to be used by the subsolve
933: - mpd - the maximum dimension allowed for the projected problem
935: Options Database Key:
936: + -eps_stoar_nev <nev> - Sets the number of eigenvalues
937: . -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
938: - -eps_stoar_mpd <mpd> - Sets the maximum projected dimension
940: Level: advanced
942: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
943: @*/
944: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
945: {
953: PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
954: return(0);
955: }
957: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
958: {
959: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
962: if (nev) *nev = ctx->nev;
963: if (ncv) *ncv = ctx->ncv;
964: if (mpd) *mpd = ctx->mpd;
965: return(0);
966: }
968: /*@
969: PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
970: step in case of doing spectrum slicing for a computational interval.
972: Not Collective
974: Input Parameter:
975: . pep - the eigenproblem solver context
977: Output Parameters:
978: + nev - number of eigenvalues to compute
979: . ncv - the maximum dimension of the subspace to be used by the subsolve
980: - mpd - the maximum dimension allowed for the projected problem
982: Level: advanced
984: .seealso: PEPSTOARSetDimensions()
985: @*/
986: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
987: {
992: PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
993: return(0);
994: }
996: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
997: {
998: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1001: ctx->checket = checket;
1002: pep->state = PEP_STATE_INITIAL;
1003: return(0);
1004: }
1006: /*@
1007: PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
1008: obtained throughout the spectrum slicing computation have the same definite type.
1010: Logically Collective on PEP
1012: Input Parameters:
1013: + pep - the eigenproblem solver context
1014: - checket - check eigenvalue type
1016: Options Database Key:
1017: . -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
1018: bool value (0/1/no/yes/true/false)
1020: Notes:
1021: This option is relevant only for spectrum slicing computations, but it is
1022: ignored if the problem type is PEP_HYPERBOLIC.
1024: This flag is turned on by default, to guarantee that the computed eigenvalues
1025: have the same type (otherwise the computed solution might be wrong). But since
1026: the check is computationally quite expensive, the check may be turned off if
1027: the user knows for sure that all eigenvalues in the requested interval have
1028: the same type.
1030: Level: advanced
1032: .seealso: PEPSetProblemType(), PEPSetInterval()
1033: @*/
1034: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
1035: {
1041: PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
1042: return(0);
1043: }
1045: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
1046: {
1047: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1050: *checket = ctx->checket;
1051: return(0);
1052: }
1054: /*@
1055: PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
1056: check in spectrum slicing.
1058: Not Collective
1060: Input Parameter:
1061: . pep - the eigenproblem solver context
1063: Output Parameter:
1064: . checket - whether eigenvalue type must be checked during spectrum slcing
1066: Level: advanced
1068: .seealso: PEPSTOARSetCheckEigenvalueType()
1069: @*/
1070: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1071: {
1077: PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1078: return(0);
1079: }
1081: PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1082: {
1084: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1085: PetscBool isascii;
1088: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1089: if (isascii) {
1090: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
1091: PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
1092: if (pep->which==PEP_ALL && !ctx->hyperbolic) {
1093: PetscViewerASCIIPrintf(viewer," checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled");
1094: }
1095: }
1096: return(0);
1097: }
1099: PetscErrorCode PEPReset_STOAR(PEP pep)
1100: {
1104: if (pep->which==PEP_ALL) {
1105: PEPReset_STOAR_QSlice(pep);
1106: }
1107: return(0);
1108: }
1110: PetscErrorCode PEPDestroy_STOAR(PEP pep)
1111: {
1113: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1116: BVDestroy(&ctx->V);
1117: PetscFree(pep->data);
1118: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL);
1119: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL);
1120: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL);
1121: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL);
1122: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL);
1123: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL);
1124: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL);
1125: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL);
1126: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL);
1127: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL);
1128: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL);
1129: return(0);
1130: }
1132: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1133: {
1135: PEP_STOAR *ctx;
1138: PetscNewLog(pep,&ctx);
1139: pep->data = (void*)ctx;
1140: ctx->lock = PETSC_TRUE;
1141: ctx->nev = 1;
1142: ctx->alpha = 1.0;
1143: ctx->beta = 0.0;
1144: ctx->checket = PETSC_TRUE;
1146: pep->ops->setup = PEPSetUp_STOAR;
1147: pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1148: pep->ops->destroy = PEPDestroy_STOAR;
1149: pep->ops->view = PEPView_STOAR;
1150: pep->ops->backtransform = PEPBackTransform_Default;
1151: pep->ops->computevectors = PEPComputeVectors_Default;
1152: pep->ops->extractvectors = PEPExtractVectors_TOAR;
1153: pep->ops->reset = PEPReset_STOAR;
1155: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR);
1156: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR);
1157: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR);
1158: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR);
1159: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR);
1160: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR);
1161: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR);
1162: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR);
1163: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR);
1164: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR);
1165: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR);
1166: return(0);
1167: }