Actual source code: stoar.c

slepc-3.13.4 2020-09-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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:   PetscErrorCode     ierr;
 52:   PEP_STOAR_MATSHELL *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: static PetscErrorCode MatDestroy_STOAR(Mat A)
 66: {
 67:   PEP_STOAR_MATSHELL *ctx;
 68:   PetscErrorCode     ierr;

 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:   PEP_STOAR_MATSHELL *ctxMat[3];
 82:   PEP_STOAR          *ctx=(PEP_STOAR*)pep->data;
 83:   PetscErrorCode     ierr;

 86:   for (i=0;i<3;i++) {
 87:     STGetMatrixTransformed(pep->st,i,&D[i]); /* D[2] = M */
 88:   }
 89:   MatGetLocalSize(D[2],&m,&n);

 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_STOAR);
 95:     MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_STOAR);
 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==PETSC_DEFAULT) 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=PETSC_FALSE,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:     /* undocumented option to use a cheaper locking instead of the true locking */
345:     PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL);
346:   }
347:   BVGetSizes(pep->V,NULL,NULL,&ld);
348:   STGetShift(pep->st,&sigma);
349:   STGetTransform(pep->st,&flg);
350:   if (pep->sfactor!=1.0) {
351:     if (!flg) {
352:       pep->target /= pep->sfactor;
353:       RGPushScale(pep->rg,1.0/pep->sfactor);
354:       STScaleShift(pep->st,1.0/pep->sfactor);
355:       sigma /= pep->sfactor;
356:     } else {
357:       PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
358:       pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
359:       RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
360:       STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
361:     }
362:   }
363:   if (flg) sigma = 0.0;

365:   /* Get the starting Arnoldi vector */
366:   BVTensorBuildFirstColumn(ctx->V,pep->nini);
367:   DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
368:   VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
369:   BVSetActiveColumns(ctx->V,0,1);
370:   BVGetSignature(ctx->V,vomega);
371:   VecGetArray(vomega,&om);
372:   omega[0] = PetscRealPart(om[0]);
373:   VecRestoreArray(vomega,&om);
374:   DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
375:   VecDestroy(&vomega);

377:   /* Restart loop */
378:   l = 0;
379:   DSGetLeadingDimension(pep->ds,&ldds);
380:   while (pep->reason == PEP_CONVERGED_ITERATING) {
381:     pep->its++;
382:     DSGetArrayReal(pep->ds,DS_MAT_T,&a);
383:     b = a+ldds;
384:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);

386:     /* Compute an nv-step Lanczos factorization */
387:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
388:     PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work);
389:     beta = b[nv-1];
390:     if (symmlost && nv==pep->nconv+l) {
391:       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
392:       pep->nconv = nconv;
393:       if (falselock || !ctx->lock) {
394:        BVSetActiveColumns(ctx->V,0,pep->nconv);
395:        BVTensorCompress(ctx->V,0);
396:       }
397:       break;
398:     }
399:     DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
400:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
401:     DSSetDimensions(pep->ds,nv,0,pep->nconv,pep->nconv+l);
402:     if (l==0) {
403:       DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
404:     } else {
405:       DSSetState(pep->ds,DS_STATE_RAW);
406:     }

408:     /* Solve projected problem */
409:     DSSolve(pep->ds,pep->eigr,pep->eigi);
410:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
411:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);

413:     /* Check convergence */
414:     /* PEPSTOARpreKConvergence(pep,nv,&norm,pep->work);*/
415:     norm = 1.0;
416:     DSGetDimensions(pep->ds,NULL,NULL,NULL,NULL,&t);
417:     PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k);
418:     (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);

420:     /* Update l */
421:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
422:     else {
423:       l = PetscMax(1,(PetscInt)((nv-k)/2));
424:       l = PetscMin(l,t);
425:       if (!breakdown) {
426:         DSGetArrayReal(pep->ds,DS_MAT_T,&a);
427:         if (*(a+ldds+k+l-1)!=0) {
428:           if (k+l<nv-1) l = l+1;
429:           else l = l-1;
430:         }
431:         /* Prepare the Rayleigh quotient for restart */
432:         DSGetArray(pep->ds,DS_MAT_Q,&Q);
433:         DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
434:         r = a + 2*ldds;
435:         for (j=k;j<k+l;j++) {
436:           r[j] = PetscRealPart(Q[nv-1+j*ldds]*beta);
437:         }
438:         b = a+ldds;
439:         b[k+l-1] = r[k+l-1];
440:         omega[k+l] = omega[nv];
441:         DSRestoreArray(pep->ds,DS_MAT_Q,&Q);
442:         DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
443:         DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
444:       }
445:     }
446:     nconv = k;
447:     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */

449:     /* Update S */
450:     DSGetMat(pep->ds,DS_MAT_Q,&MQ);
451:     BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
452:     MatDestroy(&MQ);

454:     /* Copy last column of S */
455:     BVCopyColumn(ctx->V,nv,k+l);
456:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
457:     VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
458:     VecGetArray(vomega,&om);
459:     for (j=0;j<k+l;j++) om[j] = omega[j];
460:     VecRestoreArray(vomega,&om);
461:     BVSetActiveColumns(ctx->V,0,k+l);
462:     BVSetSignature(ctx->V,vomega);
463:     VecDestroy(&vomega);
464:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);

466:     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
467:       /* stop if breakdown */
468:       PetscInfo2(pep,"Breakdown TOAR method (it=%D norm=%g)\n",pep->its,(double)beta);
469:       pep->reason = PEP_DIVERGED_BREAKDOWN;
470:     }
471:     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
472:     BVGetActiveColumns(pep->V,NULL,&nq);
473:     if (k+l+deg<=nq) {
474:       BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
475:       if (!falselock && ctx->lock) {
476:         BVTensorCompress(ctx->V,k-pep->nconv);
477:       } else {
478:         BVTensorCompress(ctx->V,0);
479:       }
480:     }
481:     pep->nconv = k;
482:     PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
483:   }

485:   if (pep->nconv>0) {
486:     BVSetActiveColumns(ctx->V,0,pep->nconv);
487:     BVGetActiveColumns(pep->V,NULL,&nq);
488:     BVSetActiveColumns(pep->V,0,nq);
489:     if (nq>pep->nconv) {
490:       BVTensorCompress(ctx->V,pep->nconv);
491:       BVSetActiveColumns(pep->V,0,pep->nconv);
492:     }
493:   }
494:   STGetTransform(pep->st,&flg);
495:   if (!flg && pep->ops->backtransform) {
496:       (*pep->ops->backtransform)(pep);
497:   }
498:   if (pep->sfactor!=1.0) {
499:     for (j=0;j<pep->nconv;j++) {
500:       pep->eigr[j] *= pep->sfactor;
501:       pep->eigi[j] *= pep->sfactor;
502:     }
503:   }
504:   /* restore original values */
505:   if (!flg) {
506:     pep->target *= pep->sfactor;
507:     STScaleShift(pep->st,pep->sfactor);
508:   } else {
509:     STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
510:     pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
511:   }
512:   if (pep->sfactor!=1.0) { RGPopScale(pep->rg); }

514:   /* truncate Schur decomposition and change the state to raw so that
515:      DSVectors() computes eigenvectors from scratch */
516:   DSSetDimensions(pep->ds,pep->nconv,0,0,0);
517:   DSSetState(pep->ds,DS_STATE_RAW);
518:   return(0);
519: }

521: PetscErrorCode PEPSetFromOptions_STOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
522: {
524:   PetscBool      flg,lock,b,f1,f2,f3;
525:   PetscInt       i,j,k;
526:   PetscReal      array[2]={0,0};
527:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;

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

532:     PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg);
533:     if (flg) { PEPSTOARSetLocking(pep,lock); }

535:     b = ctx->detect;
536:     PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg);
537:     if (flg) { PEPSTOARSetDetectZeros(pep,b); }

539:     i = 1;
540:     j = k = PETSC_DECIDE;
541:     PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1);
542:     PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2);
543:     PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3);
544:     if (f1 || f2 || f3) { PEPSTOARSetDimensions(pep,i,j,k); }

546:     k = 2;
547:     PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg);
548:     if (flg) {
549:       PEPSTOARSetLinearization(pep,array[0],array[1]);
550:     }

552:     b = ctx->checket;
553:     PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg);
554:     if (flg) { PEPSTOARSetCheckEigenvalueType(pep,b); }

556:   PetscOptionsTail();
557:   return(0);
558: }

560: static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
561: {
562:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

565:   ctx->lock = lock;
566:   return(0);
567: }

569: /*@
570:    PEPSTOARSetLocking - Choose between locking and non-locking variants of
571:    the STOAR method.

573:    Logically Collective on pep

575:    Input Parameters:
576: +  pep  - the eigenproblem solver context
577: -  lock - true if the locking variant must be selected

579:    Options Database Key:
580: .  -pep_stoar_locking - Sets the locking flag

582:    Notes:
583:    The default is to lock converged eigenpairs when the method restarts.
584:    This behaviour can be changed so that all directions are kept in the
585:    working subspace even if already converged to working accuracy (the
586:    non-locking variant).

588:    Level: advanced

590: .seealso: PEPSTOARGetLocking()
591: @*/
592: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
593: {

599:   PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
600:   return(0);
601: }

603: static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
604: {
605:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

608:   *lock = ctx->lock;
609:   return(0);
610: }

612: /*@
613:    PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.

615:    Not Collective

617:    Input Parameter:
618: .  pep - the eigenproblem solver context

620:    Output Parameter:
621: .  lock - the locking flag

623:    Level: advanced

625: .seealso: PEPSTOARSetLocking()
626: @*/
627: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
628: {

634:   PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
635:   return(0);
636: }

638: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
639: {
641:   PetscInt       i,numsh;
642:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
643:   PEP_SR         sr = ctx->sr;

646:   if (!pep->state) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
647:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
648:   switch (pep->state) {
649:   case PEP_STATE_INITIAL:
650:     break;
651:   case PEP_STATE_SETUP:
652:     if (n) *n = 2;
653:     if (shifts) {
654:       PetscMalloc1(2,shifts);
655:       (*shifts)[0] = pep->inta;
656:       (*shifts)[1] = pep->intb;
657:     }
658:     if (inertias) {
659:       PetscMalloc1(2,inertias);
660:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
661:       (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
662:     }
663:     break;
664:   case PEP_STATE_SOLVED:
665:   case PEP_STATE_EIGENVECTORS:
666:     numsh = ctx->nshifts;
667:     if (n) *n = numsh;
668:     if (shifts) {
669:       PetscMalloc1(numsh,shifts);
670:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
671:     }
672:     if (inertias) {
673:       PetscMalloc1(numsh,inertias);
674:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
675:     }
676:     break;
677:   }
678:   return(0);
679: }

681: /*@C
682:    PEPSTOARGetInertias - Gets the values of the shifts and their
683:    corresponding inertias in case of doing spectrum slicing for a
684:    computational interval.

686:    Not Collective

688:    Input Parameter:
689: .  pep - the eigenproblem solver context

691:    Output Parameters:
692: +  n        - number of shifts, including the endpoints of the interval
693: .  shifts   - the values of the shifts used internally in the solver
694: -  inertias - the values of the inertia in each shift

696:    Notes:
697:    If called after PEPSolve(), all shifts used internally by the solver are
698:    returned (including both endpoints and any intermediate ones). If called
699:    before PEPSolve() and after PEPSetUp() then only the information of the
700:    endpoints of subintervals is available.

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

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

707:    Fortran Notes:
708:    The calling sequence from Fortran is
709: .vb
710:    PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
711:    integer n
712:    double precision shifts(*)
713:    integer inertias(*)
714: .ve
715:    The arrays should be at least of length n. The value of n can be determined
716:    by an initial call
717: .vb
718:    PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
719: .ve

721:    Level: advanced

723: .seealso: PEPSetInterval()
724: @*/
725: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
726: {

732:   PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
733:   return(0);
734: }

736: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
737: {
738:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

741:   ctx->detect = detect;
742:   pep->state  = PEP_STATE_INITIAL;
743:   return(0);
744: }

746: /*@
747:    PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
748:    zeros during the factorizations throughout the spectrum slicing computation.

750:    Logically Collective on pep

752:    Input Parameters:
753: +  pep    - the eigenproblem solver context
754: -  detect - check for zeros

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

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

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

767:    Level: advanced

769: .seealso: PEPSetInterval()
770: @*/
771: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
772: {

778:   PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
779:   return(0);
780: }

782: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
783: {
784:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

787:   *detect = ctx->detect;
788:   return(0);
789: }

791: /*@
792:    PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
793:    in spectrum slicing.

795:    Not Collective

797:    Input Parameter:
798: .  pep - the eigenproblem solver context

800:    Output Parameter:
801: .  detect - whether zeros detection is enforced during factorizations

803:    Level: advanced

805: .seealso: PEPSTOARSetDetectZeros()
806: @*/
807: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
808: {

814:   PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
815:   return(0);
816: }

818: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
819: {
820:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

823:   if (beta==0.0 && alpha==0.0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
824:   ctx->alpha = alpha;
825:   ctx->beta  = beta;
826:   return(0);
827: }

829: /*@
830:    PEPSTOARSetLinearization - Set the coefficients that define
831:    the linearization of a quadratic eigenproblem.

833:    Logically Collective on pep

835:    Input Parameters:
836: +  pep   - polynomial eigenvalue solver
837: .  alpha - first parameter of the linearization
838: -  beta  - second parameter of the linearization

840:    Options Database Key:
841: .  -pep_stoar_linearization <alpha,beta> - Sets the coefficients

843:    Notes:
844:    Cannot pass zero for both alpha and beta. The default values are
845:    alpha=1 and beta=0.

847:    Level: advanced

849: .seealso: PEPSTOARGetLinearization()
850: @*/
851: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
852: {

859:   PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
860:   return(0);
861: }

863: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
864: {
865:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

868:   if (alpha) *alpha = ctx->alpha;
869:   if (beta)  *beta  = ctx->beta;
870:   return(0);
871: }

873: /*@
874:    PEPSTOARGetLinearization - Returns the coefficients that define
875:    the linearization of a quadratic eigenproblem.

877:    Not Collective

879:    Input Parameter:
880: .  pep  - polynomial eigenvalue solver

882:    Output Parameters:
883: +  alpha - the first parameter of the linearization
884: -  beta  - the second parameter of the linearization

886:    Level: advanced

888: .seealso: PEPSTOARSetLinearization()
889: @*/
890: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
891: {

896:   PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
897:   return(0);
898: }

900: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
901: {
902:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

905:   if (nev<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
906:   ctx->nev = nev;
907:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
908:     ctx->ncv = PETSC_DEFAULT;
909:   } else {
910:     if (ncv<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
911:     ctx->ncv = ncv;
912:   }
913:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
914:     ctx->mpd = PETSC_DEFAULT;
915:   } else {
916:     if (mpd<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
917:     ctx->mpd = mpd;
918:   }
919:   pep->state = PEP_STATE_INITIAL;
920:   return(0);
921: }

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

928:    Logically Collective on pep

930:    Input Parameters:
931: +  pep - the eigenproblem solver context
932: .  nev - number of eigenvalues to compute
933: .  ncv - the maximum dimension of the subspace to be used by the subsolve
934: -  mpd - the maximum dimension allowed for the projected problem

936:    Options Database Key:
937: +  -eps_stoar_nev <nev> - Sets the number of eigenvalues
938: .  -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
939: -  -eps_stoar_mpd <mpd> - Sets the maximum projected dimension

941:    Level: advanced

943: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
944: @*/
945: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
946: {

954:   PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
955:   return(0);
956: }

958: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
959: {
960:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

963:   if (nev) *nev = ctx->nev;
964:   if (ncv) *ncv = ctx->ncv;
965:   if (mpd) *mpd = ctx->mpd;
966:   return(0);
967: }

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

973:    Not Collective

975:    Input Parameter:
976: .  pep - the eigenproblem solver context

978:    Output Parameters:
979: +  nev - number of eigenvalues to compute
980: .  ncv - the maximum dimension of the subspace to be used by the subsolve
981: -  mpd - the maximum dimension allowed for the projected problem

983:    Level: advanced

985: .seealso: PEPSTOARSetDimensions()
986: @*/
987: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
988: {

993:   PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
994:   return(0);
995: }

997: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
998: {
999:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

1002:   ctx->checket = checket;
1003:   pep->state   = PEP_STATE_INITIAL;
1004:   return(0);
1005: }

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

1011:    Logically Collective on pep

1013:    Input Parameters:
1014: +  pep     - the eigenproblem solver context
1015: -  checket - check eigenvalue type

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

1021:    Notes:
1022:    This option is relevant only for spectrum slicing computations, but it is
1023:    ignored if the problem type is PEP_HYPERBOLIC.

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

1031:    Level: advanced

1033: .seealso: PEPSetProblemType(), PEPSetInterval()
1034: @*/
1035: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
1036: {

1042:   PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
1043:   return(0);
1044: }

1046: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
1047: {
1048:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

1051:   *checket = ctx->checket;
1052:   return(0);
1053: }

1055: /*@
1056:    PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
1057:    check in spectrum slicing.

1059:    Not Collective

1061:    Input Parameter:
1062: .  pep - the eigenproblem solver context

1064:    Output Parameter:
1065: .  checket - whether eigenvalue type must be checked during spectrum slcing

1067:    Level: advanced

1069: .seealso: PEPSTOARSetCheckEigenvalueType()
1070: @*/
1071: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1072: {

1078:   PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1079:   return(0);
1080: }

1082: PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1083: {
1085:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
1086:   PetscBool      isascii;

1089:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1090:   if (isascii) {
1091:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1092:     PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
1093:     if (pep->which==PEP_ALL && !ctx->hyperbolic) {
1094:       PetscViewerASCIIPrintf(viewer,"  checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled");
1095:     }
1096:   }
1097:   return(0);
1098: }

1100: PetscErrorCode PEPReset_STOAR(PEP pep)
1101: {

1105:   if (pep->which==PEP_ALL) {
1106:     PEPReset_STOAR_QSlice(pep);
1107:   }
1108:   return(0);
1109: }

1111: PetscErrorCode PEPDestroy_STOAR(PEP pep)
1112: {
1114:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;

1117:   BVDestroy(&ctx->V);
1118:   PetscFree(pep->data);
1119:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL);
1120:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL);
1121:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL);
1122:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL);
1123:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL);
1124:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL);
1125:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL);
1126:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL);
1127:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL);
1128:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL);
1129:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL);
1130:   return(0);
1131: }

1133: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1134: {
1136:   PEP_STOAR      *ctx;

1139:   PetscNewLog(pep,&ctx);
1140:   pep->data = (void*)ctx;
1141:   ctx->lock    = PETSC_TRUE;
1142:   ctx->nev     = 1;
1143:   ctx->ncv     = PETSC_DEFAULT;
1144:   ctx->mpd     = PETSC_DEFAULT;
1145:   ctx->alpha   = 1.0;
1146:   ctx->beta    = 0.0;
1147:   ctx->checket = PETSC_TRUE;

1149:   pep->ops->setup          = PEPSetUp_STOAR;
1150:   pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1151:   pep->ops->destroy        = PEPDestroy_STOAR;
1152:   pep->ops->view           = PEPView_STOAR;
1153:   pep->ops->backtransform  = PEPBackTransform_Default;
1154:   pep->ops->computevectors = PEPComputeVectors_Default;
1155:   pep->ops->extractvectors = PEPExtractVectors_TOAR;
1156:   pep->ops->reset          = PEPReset_STOAR;

1158:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR);
1159:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR);
1160:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR);
1161:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR);
1162:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR);
1163:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR);
1164:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR);
1165:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR);
1166:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR);
1167:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR);
1168:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR);
1169:   return(0);
1170: }