Actual source code: stoar.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  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 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>
 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: } PEP_STOAR_MATSHELL;

 49: static PetscErrorCode MatMult_STOAR(Mat A,Vec x,Vec y)
 50: {
 51:   PEP_STOAR_MATSHELL *ctx;

 53:   MatShellGetContext(A,&ctx);
 54:   MatMult(ctx->A[0],x,y);
 55:   VecScale(y,ctx->scal[0]);
 56:   if (ctx->scal[1]) {
 57:     MatMult(ctx->A[1],x,ctx->t);
 58:     VecAXPY(y,ctx->scal[1],ctx->t);
 59:   }
 60:   PetscFunctionReturn(0);
 61: }

 63: static PetscErrorCode MatDestroy_STOAR(Mat A)
 64: {
 65:   PEP_STOAR_MATSHELL *ctx;

 67:   MatShellGetContext(A,&ctx);
 68:   VecDestroy(&ctx->t);
 69:   PetscFree(ctx);
 70:   PetscFunctionReturn(0);
 71: }

 73: PetscErrorCode PEPSTOARSetUpInnerMatrix(PEP pep,Mat *B)
 74: {
 75:   Mat                pB[4],Bs[3],D[3];
 76:   PetscInt           i,j,n,m;
 77:   PEP_STOAR_MATSHELL *ctxMat[3];
 78:   PEP_STOAR          *ctx=(PEP_STOAR*)pep->data;

 80:   for (i=0;i<3;i++) {
 81:     STGetMatrixTransformed(pep->st,i,&D[i]); /* D[2] = M */
 82:   }
 83:   MatGetLocalSize(D[2],&m,&n);

 85:   for (j=0;j<3;j++) {
 86:     PetscNew(ctxMat+j);
 87:     MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&Bs[j]);
 88:     MatShellSetOperation(Bs[j],MATOP_MULT,(void(*)(void))MatMult_STOAR);
 89:     MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_STOAR);
 90:   }
 91:   for (i=0;i<4;i++) pB[i] = NULL;
 92:   if (ctx->alpha) {
 93:     ctxMat[0]->A[0] = D[0]; ctxMat[0]->scal[0] = ctx->alpha; ctxMat[0]->scal[1] = 0.0;
 94:     ctxMat[2]->A[0] = D[2]; ctxMat[2]->scal[0] = -ctx->alpha*pep->sfactor*pep->sfactor; ctxMat[2]->scal[1] = 0.0;
 95:     pB[0] = Bs[0]; pB[3] = Bs[2];
 96:   }
 97:   if (ctx->beta) {
 98:     i = (ctx->alpha)?1:0;
 99:     ctxMat[0]->scal[1] = 0.0;
100:     ctxMat[0]->A[i] = D[1]; ctxMat[0]->scal[i] = -ctx->beta*pep->sfactor;
101:     ctxMat[1]->A[0] = D[2]; ctxMat[1]->scal[0] = -ctx->beta*pep->sfactor*pep->sfactor; ctxMat[1]->scal[1] = 0.0;
102:     pB[0] = Bs[0]; pB[1] = pB[2] = Bs[1];
103:   }
104:   BVCreateVec(pep->V,&ctxMat[0]->t);
105:   MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pB,B);
106:   for (j=0;j<3;j++) MatDestroy(&Bs[j]);
107:   PetscFunctionReturn(0);
108: }

110: PetscErrorCode PEPSetUp_STOAR(PEP pep)
111: {
112:   PetscBool         sinv,flg;
113:   PEP_STOAR         *ctx = (PEP_STOAR*)pep->data;
114:   PetscInt          ld,i;
115:   PetscReal         eta;
116:   BVOrthogType      otype;
117:   BVOrthogBlockType obtype;

119:   PEPCheckHermitian(pep);
120:   PEPCheckQuadratic(pep);
121:   PEPCheckShiftSinvert(pep);
122:   /* spectrum slicing requires special treatment of default values */
123:   if (pep->which==PEP_ALL) {
124:     pep->ops->solve = PEPSolve_STOAR_QSlice;
125:     pep->ops->extractvectors = NULL;
126:     pep->ops->setdefaultst   = NULL;
127:     PEPSetUp_STOAR_QSlice(pep);
128:   } else {
129:     PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
131:     if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
132:     pep->ops->solve = PEPSolve_STOAR;
133:     ld   = pep->ncv+2;
134:     DSSetType(pep->ds,DSGHIEP);
135:     DSSetCompact(pep->ds,PETSC_TRUE);
136:     DSSetExtraRow(pep->ds,PETSC_TRUE);
137:     DSAllocate(pep->ds,ld);
138:     PEPBasisCoefficients(pep,pep->pbc);
139:     STGetTransform(pep->st,&flg);
140:     if (!flg) {
141:       PetscFree(pep->solvematcoeffs);
142:       PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
143:       PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
144:       PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
145:       if (sinv) PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
146:       else {
147:         for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
148:         pep->solvematcoeffs[pep->nmat-1] = 1.0;
149:       }
150:     }
151:   }
152:   if (!pep->which) PEPSetWhichEigenpairs_Default(pep);
153:   PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_REGION);

155:   PEPAllocateSolution(pep,2);
156:   PEPSetWorkVecs(pep,4);
157:   BVDestroy(&ctx->V);
158:   BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
159:   BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype);
160:   BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
161:   PetscFunctionReturn(0);
162: }

164: /*
165:   Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
166: */
167: static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
168: {
169:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
170:   PetscInt       i,j,m=*M,l,lock;
171:   PetscInt       lds,d,ld,offq,nqt,ldds;
172:   Vec            v=t_[0],t=t_[1],q=t_[2];
173:   PetscReal      norm,sym=0.0,fro=0.0,*f;
174:   PetscScalar    *y,*S,*x,sigma;
175:   PetscBLASInt   j_,one=1;
176:   PetscBool      lindep,flg,sinvert=PETSC_FALSE;
177:   Mat            MS;

179:   PetscMalloc1(*M,&y);
180:   BVGetSizes(pep->V,NULL,NULL,&ld);
181:   BVTensorGetDegree(ctx->V,&d);
182:   BVGetActiveColumns(pep->V,&lock,&nqt);
183:   lds = d*ld;
184:   offq = ld;
185:   DSGetLeadingDimension(pep->ds,&ldds);
186:   *breakdown = PETSC_FALSE; /* ----- */
187:   DSGetDimensions(pep->ds,NULL,&l,NULL,NULL);
188:   BVSetActiveColumns(ctx->V,0,m);
189:   BVSetActiveColumns(pep->V,0,nqt);
190:   STGetTransform(pep->st,&flg);
191:   if (!flg) {
192:     /* spectral transformation handled by the solver */
193:     PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,"");
195:     PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert);
196:     STGetShift(pep->st,&sigma);
197:   }
198:   for (j=k;j<m;j++) {
199:     /* apply operator */
200:     BVTensorGetFactors(ctx->V,NULL,&MS);
201:     MatDenseGetArray(MS,&S);
202:     BVGetColumn(pep->V,nqt,&t);
203:     BVMultVec(pep->V,1.0,0.0,v,S+j*lds);
204:     if (!sinvert) {
205:       STMatMult(pep->st,0,v,q);
206:       BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
207:       STMatMult(pep->st,1,v,t);
208:       VecAXPY(q,pep->sfactor,t);
209:       if (ctx->beta && ctx->alpha) {
210:         STMatMult(pep->st,2,v,t);
211:         VecAXPY(q,-pep->sfactor*pep->sfactor*ctx->beta/ctx->alpha,t);
212:       }
213:       STMatSolve(pep->st,q,t);
214:       VecScale(t,-1.0/(pep->sfactor*pep->sfactor));
215:     } else {
216:       STMatMult(pep->st,1,v,q);
217:       STMatMult(pep->st,2,v,t);
218:       VecAXPY(q,sigma*pep->sfactor,t);
219:       VecScale(q,pep->sfactor);
220:       BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
221:       STMatMult(pep->st,2,v,t);
222:       VecAXPY(q,pep->sfactor*pep->sfactor,t);
223:       STMatSolve(pep->st,q,t);
224:       VecScale(t,-1.0);
225:     }
226:     BVRestoreColumn(pep->V,nqt,&t);

228:     /* orthogonalize */
229:     if (!sinvert) x = S+offq+(j+1)*lds;
230:     else x = S+(j+1)*lds;
231:     BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep);

233:     if (!lindep) {
234:       if (!sinvert) *(S+offq+(j+1)*lds+nqt) = norm;
235:       else *(S+(j+1)*lds+nqt) = norm;
236:       BVScaleColumn(pep->V,nqt,1.0/norm);
237:       nqt++;
238:     }
239:     if (!sinvert) {
240:       for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
241:       if (ctx->beta && ctx->alpha) {
242:         for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+offq+i) -= *(S+(j+1)*lds+i)*ctx->beta/ctx->alpha;
243:       }
244:     } else for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
245:     BVSetActiveColumns(pep->V,0,nqt);
246:     MatDenseRestoreArray(MS,&S);
247:     BVTensorRestoreFactors(ctx->V,NULL,&MS);

249:     /* level-2 orthogonalization */
250:     BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep);
251:     a[j] = PetscRealPart(y[j]);
252:     omega[j+1] = (norm > 0)?1.0:-1.0;
253:     BVScaleColumn(ctx->V,j+1,1.0/norm);
254:     b[j] = PetscAbsReal(norm);

256:     /* check symmetry */
257:     DSGetArrayReal(pep->ds,DS_MAT_T,&f);
258:     if (j==k) {
259:       for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
260:       for (i=0;i<l;i++) y[i] = 0.0;
261:     }
262:     DSRestoreArrayReal(pep->ds,DS_MAT_T,&f);
263:     if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
264:     PetscBLASIntCast(j,&j_);
265:     sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
266:     fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
267:     if (j>0) fro = SlepcAbs(fro,b[j-1]);
268:     if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
269:       *symmlost = PETSC_TRUE;
270:       *M=j;
271:       break;
272:     }
273:   }
274:   BVSetActiveColumns(pep->V,lock,nqt);
275:   BVSetActiveColumns(ctx->V,0,*M);
276:   PetscFree(y);
277:   PetscFunctionReturn(0);
278: }

280: #if 0
281: static PetscErrorCode PEPSTOARpreKConvergence(PEP pep,PetscInt nv,PetscReal *norm,Vec *w)
282: {
283:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
284:   PetscBLASInt   n_,one=1;
285:   PetscInt       lds=2*ctx->ld;
286:   PetscReal      t1,t2;
287:   PetscScalar    *S=ctx->S;

289:   PetscBLASIntCast(nv+2,&n_);
290:   t1 = BLASnrm2_(&n_,S+nv*2*ctx->ld,&one);
291:   t2 = BLASnrm2_(&n_,S+(nv*2+1)*ctx->ld,&one);
292:   *norm = SlepcAbs(t1,t2);
293:   BVSetActiveColumns(pep->V,0,nv+2);
294:   BVMultVec(pep->V,1.0,0.0,w[1],S+nv*lds);
295:   STMatMult(pep->st,0,w[1],w[2]);
296:   VecNorm(w[2],NORM_2,&t1);
297:   BVMultVec(pep->V,1.0,0.0,w[1],S+ctx->ld+nv*lds);
298:   STMatMult(pep->st,2,w[1],w[2]);
299:   VecNorm(w[2],NORM_2,&t2);
300:   t2 *= pep->sfactor*pep->sfactor;
301:   *norm = PetscMax(*norm,SlepcAbs(t1,t2));
302:   PetscFunctionReturn(0);
303: }
304: #endif

306: PetscErrorCode PEPSolve_STOAR(PEP pep)
307: {
308:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
309:   PetscInt       j,k,l,nv=0,ld,ldds,t,nq=0;
310:   PetscInt       nconv=0,deg=pep->nmat-1;
311:   PetscScalar    *om,sigma;
312:   PetscReal      beta,norm=1.0,*omega,*a,*b;
313:   PetscBool      breakdown,symmlost=PETSC_FALSE,sinv=PETSC_FALSE,falselock=PETSC_TRUE,flg;
314:   Mat            MQ,A;
315:   Vec            vomega;

317:   PetscCitationsRegister(citation,&cited);
318:   PEPSTOARSetUpInnerMatrix(pep,&A);
319:   BVSetMatrix(ctx->V,A,PETSC_TRUE);
320:   MatDestroy(&A);
321:   if (ctx->lock) {
322:     /* undocumented option to use a cheaper locking instead of the true locking */
323:     PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL);
324:   }
325:   BVGetSizes(pep->V,NULL,NULL,&ld);
326:   STGetShift(pep->st,&sigma);
327:   STGetTransform(pep->st,&flg);
328:   if (pep->sfactor!=1.0) {
329:     if (!flg) {
330:       pep->target /= pep->sfactor;
331:       RGPushScale(pep->rg,1.0/pep->sfactor);
332:       STScaleShift(pep->st,1.0/pep->sfactor);
333:       sigma /= pep->sfactor;
334:     } else {
335:       PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
336:       pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
337:       RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
338:       STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
339:     }
340:   }
341:   if (flg) sigma = 0.0;

343:   /* Get the starting Arnoldi vector */
344:   BVTensorBuildFirstColumn(ctx->V,pep->nini);
345:   DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
346:   VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
347:   BVSetActiveColumns(ctx->V,0,1);
348:   BVGetSignature(ctx->V,vomega);
349:   VecGetArray(vomega,&om);
350:   omega[0] = PetscRealPart(om[0]);
351:   VecRestoreArray(vomega,&om);
352:   DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
353:   VecDestroy(&vomega);

355:   /* Restart loop */
356:   l = 0;
357:   DSGetLeadingDimension(pep->ds,&ldds);
358:   while (pep->reason == PEP_CONVERGED_ITERATING) {
359:     pep->its++;
360:     DSGetArrayReal(pep->ds,DS_MAT_T,&a);
361:     b = a+ldds;
362:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);

364:     /* Compute an nv-step Lanczos factorization */
365:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
366:     PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work);
367:     beta = b[nv-1];
368:     if (symmlost && nv==pep->nconv+l) {
369:       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
370:       pep->nconv = nconv;
371:       if (falselock || !ctx->lock) {
372:        BVSetActiveColumns(ctx->V,0,pep->nconv);
373:        BVTensorCompress(ctx->V,0);
374:       }
375:       break;
376:     }
377:     DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
378:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
379:     DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l);
380:     if (l==0) DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
381:     else DSSetState(pep->ds,DS_STATE_RAW);

383:     /* Solve projected problem */
384:     DSSolve(pep->ds,pep->eigr,pep->eigi);
385:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
386:     DSUpdateExtraRow(pep->ds);
387:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);

389:     /* Check convergence */
390:     /* PEPSTOARpreKConvergence(pep,nv,&norm,pep->work);*/
391:     norm = 1.0;
392:     DSGetDimensions(pep->ds,NULL,NULL,NULL,&t);
393:     PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k);
394:     (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);

396:     /* Update l */
397:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
398:     else {
399:       l = PetscMax(1,(PetscInt)((nv-k)/2));
400:       l = PetscMin(l,t);
401:       DSGetTruncateSize(pep->ds,k,t,&l);
402:       if (!breakdown) {
403:         /* Prepare the Rayleigh quotient for restart */
404:         DSTruncate(pep->ds,k+l,PETSC_FALSE);
405:       }
406:     }
407:     nconv = k;
408:     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
409:     if (l) PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

411:     /* Update S */
412:     DSGetMat(pep->ds,DS_MAT_Q,&MQ);
413:     BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
414:     MatDestroy(&MQ);

416:     /* Copy last column of S */
417:     BVCopyColumn(ctx->V,nv,k+l);
418:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
419:     VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
420:     VecGetArray(vomega,&om);
421:     for (j=0;j<k+l;j++) om[j] = omega[j];
422:     VecRestoreArray(vomega,&om);
423:     BVSetActiveColumns(ctx->V,0,k+l);
424:     BVSetSignature(ctx->V,vomega);
425:     VecDestroy(&vomega);
426:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);

428:     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
429:       /* stop if breakdown */
430:       PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta);
431:       pep->reason = PEP_DIVERGED_BREAKDOWN;
432:     }
433:     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
434:     BVGetActiveColumns(pep->V,NULL,&nq);
435:     if (k+l+deg<=nq) {
436:       BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
437:       if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-pep->nconv);
438:       else BVTensorCompress(ctx->V,0);
439:     }
440:     pep->nconv = k;
441:     PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
442:   }

444:   if (pep->nconv>0) {
445:     BVSetActiveColumns(ctx->V,0,pep->nconv);
446:     BVGetActiveColumns(pep->V,NULL,&nq);
447:     BVSetActiveColumns(pep->V,0,nq);
448:     if (nq>pep->nconv) {
449:       BVTensorCompress(ctx->V,pep->nconv);
450:       BVSetActiveColumns(pep->V,0,pep->nconv);
451:     }
452:   }
453:   STGetTransform(pep->st,&flg);
454:   if (!flg && pep->ops->backtransform) (*pep->ops->backtransform)(pep);
455:   if (pep->sfactor!=1.0) {
456:     for (j=0;j<pep->nconv;j++) {
457:       pep->eigr[j] *= pep->sfactor;
458:       pep->eigi[j] *= pep->sfactor;
459:     }
460:   }
461:   /* restore original values */
462:   if (!flg) {
463:     pep->target *= pep->sfactor;
464:     STScaleShift(pep->st,pep->sfactor);
465:   } else {
466:     STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
467:     pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
468:   }
469:   if (pep->sfactor!=1.0) RGPopScale(pep->rg);

471:   DSTruncate(pep->ds,pep->nconv,PETSC_TRUE);
472:   PetscFunctionReturn(0);
473: }

475: PetscErrorCode PEPSetFromOptions_STOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
476: {
477:   PetscBool      flg,lock,b,f1,f2,f3;
478:   PetscInt       i,j,k;
479:   PetscReal      array[2]={0,0};
480:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;

482:   PetscOptionsHead(PetscOptionsObject,"PEP STOAR Options");

484:     PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg);
485:     if (flg) PEPSTOARSetLocking(pep,lock);

487:     b = ctx->detect;
488:     PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg);
489:     if (flg) PEPSTOARSetDetectZeros(pep,b);

491:     i = 1;
492:     j = k = PETSC_DECIDE;
493:     PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1);
494:     PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2);
495:     PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3);
496:     if (f1 || f2 || f3) PEPSTOARSetDimensions(pep,i,j,k);

498:     k = 2;
499:     PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg);
500:     if (flg) PEPSTOARSetLinearization(pep,array[0],array[1]);

502:     b = ctx->checket;
503:     PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg);
504:     if (flg) PEPSTOARSetCheckEigenvalueType(pep,b);

506:   PetscOptionsTail();
507:   PetscFunctionReturn(0);
508: }

510: static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
511: {
512:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

514:   ctx->lock = lock;
515:   PetscFunctionReturn(0);
516: }

518: /*@
519:    PEPSTOARSetLocking - Choose between locking and non-locking variants of
520:    the STOAR method.

522:    Logically Collective on pep

524:    Input Parameters:
525: +  pep  - the eigenproblem solver context
526: -  lock - true if the locking variant must be selected

528:    Options Database Key:
529: .  -pep_stoar_locking - Sets the locking flag

531:    Notes:
532:    The default is to lock converged eigenpairs when the method restarts.
533:    This behaviour can be changed so that all directions are kept in the
534:    working subspace even if already converged to working accuracy (the
535:    non-locking variant).

537:    Level: advanced

539: .seealso: PEPSTOARGetLocking()
540: @*/
541: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
542: {
545:   PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
546:   PetscFunctionReturn(0);
547: }

549: static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
550: {
551:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

553:   *lock = ctx->lock;
554:   PetscFunctionReturn(0);
555: }

557: /*@
558:    PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.

560:    Not Collective

562:    Input Parameter:
563: .  pep - the eigenproblem solver context

565:    Output Parameter:
566: .  lock - the locking flag

568:    Level: advanced

570: .seealso: PEPSTOARSetLocking()
571: @*/
572: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
573: {
576:   PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
577:   PetscFunctionReturn(0);
578: }

580: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
581: {
582:   PetscInt       i,numsh;
583:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
584:   PEP_SR         sr = ctx->sr;

588:   switch (pep->state) {
589:   case PEP_STATE_INITIAL:
590:     break;
591:   case PEP_STATE_SETUP:
592:     if (n) *n = 2;
593:     if (shifts) {
594:       PetscMalloc1(2,shifts);
595:       (*shifts)[0] = pep->inta;
596:       (*shifts)[1] = pep->intb;
597:     }
598:     if (inertias) {
599:       PetscMalloc1(2,inertias);
600:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
601:       (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
602:     }
603:     break;
604:   case PEP_STATE_SOLVED:
605:   case PEP_STATE_EIGENVECTORS:
606:     numsh = ctx->nshifts;
607:     if (n) *n = numsh;
608:     if (shifts) {
609:       PetscMalloc1(numsh,shifts);
610:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
611:     }
612:     if (inertias) {
613:       PetscMalloc1(numsh,inertias);
614:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
615:     }
616:     break;
617:   }
618:   PetscFunctionReturn(0);
619: }

621: /*@C
622:    PEPSTOARGetInertias - Gets the values of the shifts and their
623:    corresponding inertias in case of doing spectrum slicing for a
624:    computational interval.

626:    Not Collective

628:    Input Parameter:
629: .  pep - the eigenproblem solver context

631:    Output Parameters:
632: +  n        - number of shifts, including the endpoints of the interval
633: .  shifts   - the values of the shifts used internally in the solver
634: -  inertias - the values of the inertia in each shift

636:    Notes:
637:    If called after PEPSolve(), all shifts used internally by the solver are
638:    returned (including both endpoints and any intermediate ones). If called
639:    before PEPSolve() and after PEPSetUp() then only the information of the
640:    endpoints of subintervals is available.

642:    This function is only available for spectrum slicing runs.

644:    The returned arrays should be freed by the user. Can pass NULL in any of
645:    the two arrays if not required.

647:    Fortran Notes:
648:    The calling sequence from Fortran is
649: .vb
650:    PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
651:    integer n
652:    double precision shifts(*)
653:    integer inertias(*)
654: .ve
655:    The arrays should be at least of length n. The value of n can be determined
656:    by an initial call
657: .vb
658:    PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
659: .ve

661:    Level: advanced

663: .seealso: PEPSetInterval()
664: @*/
665: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
666: {
669:   PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
670:   PetscFunctionReturn(0);
671: }

673: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
674: {
675:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

677:   ctx->detect = detect;
678:   pep->state  = PEP_STATE_INITIAL;
679:   PetscFunctionReturn(0);
680: }

682: /*@
683:    PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
684:    zeros during the factorizations throughout the spectrum slicing computation.

686:    Logically Collective on pep

688:    Input Parameters:
689: +  pep    - the eigenproblem solver context
690: -  detect - check for zeros

692:    Options Database Key:
693: .  -pep_stoar_detect_zeros - Check for zeros; this takes an optional
694:    bool value (0/1/no/yes/true/false)

696:    Notes:
697:    A zero in the factorization indicates that a shift coincides with an eigenvalue.

699:    This flag is turned off by default, and may be necessary in some cases.
700:    This feature currently requires an external package for factorizations
701:    with support for zero detection, e.g. MUMPS.

703:    Level: advanced

705: .seealso: PEPSetInterval()
706: @*/
707: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
708: {
711:   PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
712:   PetscFunctionReturn(0);
713: }

715: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
716: {
717:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

719:   *detect = ctx->detect;
720:   PetscFunctionReturn(0);
721: }

723: /*@
724:    PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
725:    in spectrum slicing.

727:    Not Collective

729:    Input Parameter:
730: .  pep - the eigenproblem solver context

732:    Output Parameter:
733: .  detect - whether zeros detection is enforced during factorizations

735:    Level: advanced

737: .seealso: PEPSTOARSetDetectZeros()
738: @*/
739: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
740: {
743:   PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
744:   PetscFunctionReturn(0);
745: }

747: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
748: {
749:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

752:   ctx->alpha = alpha;
753:   ctx->beta  = beta;
754:   PetscFunctionReturn(0);
755: }

757: /*@
758:    PEPSTOARSetLinearization - Set the coefficients that define
759:    the linearization of a quadratic eigenproblem.

761:    Logically Collective on pep

763:    Input Parameters:
764: +  pep   - polynomial eigenvalue solver
765: .  alpha - first parameter of the linearization
766: -  beta  - second parameter of the linearization

768:    Options Database Key:
769: .  -pep_stoar_linearization <alpha,beta> - Sets the coefficients

771:    Notes:
772:    Cannot pass zero for both alpha and beta. The default values are
773:    alpha=1 and beta=0.

775:    Level: advanced

777: .seealso: PEPSTOARGetLinearization()
778: @*/
779: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
780: {
784:   PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
785:   PetscFunctionReturn(0);
786: }

788: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
789: {
790:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

792:   if (alpha) *alpha = ctx->alpha;
793:   if (beta)  *beta  = ctx->beta;
794:   PetscFunctionReturn(0);
795: }

797: /*@
798:    PEPSTOARGetLinearization - Returns the coefficients that define
799:    the linearization of a quadratic eigenproblem.

801:    Not Collective

803:    Input Parameter:
804: .  pep  - polynomial eigenvalue solver

806:    Output Parameters:
807: +  alpha - the first parameter of the linearization
808: -  beta  - the second parameter of the linearization

810:    Level: advanced

812: .seealso: PEPSTOARSetLinearization()
813: @*/
814: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
815: {
817:   PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
818:   PetscFunctionReturn(0);
819: }

821: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
822: {
823:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

826:   ctx->nev = nev;
827:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
828:     ctx->ncv = PETSC_DEFAULT;
829:   } else {
831:     ctx->ncv = ncv;
832:   }
833:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
834:     ctx->mpd = PETSC_DEFAULT;
835:   } else {
837:     ctx->mpd = mpd;
838:   }
839:   pep->state = PEP_STATE_INITIAL;
840:   PetscFunctionReturn(0);
841: }

843: /*@
844:    PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
845:    step in case of doing spectrum slicing for a computational interval.
846:    The meaning of the parameters is the same as in PEPSetDimensions().

848:    Logically Collective on pep

850:    Input Parameters:
851: +  pep - the eigenproblem solver context
852: .  nev - number of eigenvalues to compute
853: .  ncv - the maximum dimension of the subspace to be used by the subsolve
854: -  mpd - the maximum dimension allowed for the projected problem

856:    Options Database Key:
857: +  -eps_stoar_nev <nev> - Sets the number of eigenvalues
858: .  -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
859: -  -eps_stoar_mpd <mpd> - Sets the maximum projected dimension

861:    Level: advanced

863: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
864: @*/
865: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
866: {
871:   PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
872:   PetscFunctionReturn(0);
873: }

875: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
876: {
877:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

879:   if (nev) *nev = ctx->nev;
880:   if (ncv) *ncv = ctx->ncv;
881:   if (mpd) *mpd = ctx->mpd;
882:   PetscFunctionReturn(0);
883: }

885: /*@
886:    PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
887:    step in case of doing spectrum slicing for a computational interval.

889:    Not Collective

891:    Input Parameter:
892: .  pep - the eigenproblem solver context

894:    Output Parameters:
895: +  nev - number of eigenvalues to compute
896: .  ncv - the maximum dimension of the subspace to be used by the subsolve
897: -  mpd - the maximum dimension allowed for the projected problem

899:    Level: advanced

901: .seealso: PEPSTOARSetDimensions()
902: @*/
903: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
904: {
906:   PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
907:   PetscFunctionReturn(0);
908: }

910: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
911: {
912:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

914:   ctx->checket = checket;
915:   pep->state   = PEP_STATE_INITIAL;
916:   PetscFunctionReturn(0);
917: }

919: /*@
920:    PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
921:    obtained throughout the spectrum slicing computation have the same definite type.

923:    Logically Collective on pep

925:    Input Parameters:
926: +  pep     - the eigenproblem solver context
927: -  checket - check eigenvalue type

929:    Options Database Key:
930: .  -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
931:    bool value (0/1/no/yes/true/false)

933:    Notes:
934:    This option is relevant only for spectrum slicing computations, but it is
935:    ignored if the problem type is PEP_HYPERBOLIC.

937:    This flag is turned on by default, to guarantee that the computed eigenvalues
938:    have the same type (otherwise the computed solution might be wrong). But since
939:    the check is computationally quite expensive, the check may be turned off if
940:    the user knows for sure that all eigenvalues in the requested interval have
941:    the same type.

943:    Level: advanced

945: .seealso: PEPSetProblemType(), PEPSetInterval()
946: @*/
947: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
948: {
951:   PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
952:   PetscFunctionReturn(0);
953: }

955: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
956: {
957:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

959:   *checket = ctx->checket;
960:   PetscFunctionReturn(0);
961: }

963: /*@
964:    PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
965:    check in spectrum slicing.

967:    Not Collective

969:    Input Parameter:
970: .  pep - the eigenproblem solver context

972:    Output Parameter:
973: .  checket - whether eigenvalue type must be checked during spectrum slcing

975:    Level: advanced

977: .seealso: PEPSTOARSetCheckEigenvalueType()
978: @*/
979: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
980: {
983:   PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
984:   PetscFunctionReturn(0);
985: }

987: PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
988: {
989:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
990:   PetscBool      isascii;

992:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
993:   if (isascii) {
994:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
995:     PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
996:     if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscViewerASCIIPrintf(viewer,"  checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled");
997:   }
998:   PetscFunctionReturn(0);
999: }

1001: PetscErrorCode PEPReset_STOAR(PEP pep)
1002: {
1003:   if (pep->which==PEP_ALL) PEPReset_STOAR_QSlice(pep);
1004:   PetscFunctionReturn(0);
1005: }

1007: PetscErrorCode PEPDestroy_STOAR(PEP pep)
1008: {
1009:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;

1011:   BVDestroy(&ctx->V);
1012:   PetscFree(pep->data);
1013:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL);
1014:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL);
1015:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL);
1016:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL);
1017:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL);
1018:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL);
1019:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL);
1020:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL);
1021:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL);
1022:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL);
1023:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL);
1024:   PetscFunctionReturn(0);
1025: }

1027: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1028: {
1029:   PEP_STOAR      *ctx;

1031:   PetscNewLog(pep,&ctx);
1032:   pep->data = (void*)ctx;

1034:   pep->lineariz = PETSC_TRUE;
1035:   ctx->lock     = PETSC_TRUE;
1036:   ctx->nev      = 1;
1037:   ctx->ncv      = PETSC_DEFAULT;
1038:   ctx->mpd      = PETSC_DEFAULT;
1039:   ctx->alpha    = 1.0;
1040:   ctx->beta     = 0.0;
1041:   ctx->checket  = PETSC_TRUE;

1043:   pep->ops->setup          = PEPSetUp_STOAR;
1044:   pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1045:   pep->ops->destroy        = PEPDestroy_STOAR;
1046:   pep->ops->view           = PEPView_STOAR;
1047:   pep->ops->backtransform  = PEPBackTransform_Default;
1048:   pep->ops->computevectors = PEPComputeVectors_Default;
1049:   pep->ops->extractvectors = PEPExtractVectors_TOAR;
1050:   pep->ops->reset          = PEPReset_STOAR;

1052:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR);
1053:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR);
1054:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR);
1055:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR);
1056:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR);
1057:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR);
1058:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR);
1059:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR);
1060:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR);
1061:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR);
1062:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR);
1063:   PetscFunctionReturn(0);
1064: }