Actual source code: stoar.c

slepc-3.22.2 2024-12-02
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:   PetscFunctionBegin;
 54:   PetscCall(MatShellGetContext(A,&ctx));
 55:   PetscCall(MatMult(ctx->A[0],x,y));
 56:   PetscCall(VecScale(y,ctx->scal[0]));
 57:   if (ctx->scal[1]) {
 58:     PetscCall(MatMult(ctx->A[1],x,ctx->t));
 59:     PetscCall(VecAXPY(y,ctx->scal[1],ctx->t));
 60:   }
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

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

 68:   PetscFunctionBegin;
 69:   PetscCall(MatShellGetContext(A,&ctx));
 70:   PetscCall(VecDestroy(&ctx->t));
 71:   PetscCall(PetscFree(ctx));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

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

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

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

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

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

158:   PetscCall(PEPAllocateSolution(pep,2));
159:   PetscCall(PEPSetWorkVecs(pep,4));
160:   PetscCall(BVDestroy(&ctx->V));
161:   PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
162:   PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
163:   PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

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

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

232:     /* orthogonalize */
233:     if (!sinvert) x = S+offq+(j+1)*lds;
234:     else x = S+(j+1)*lds;
235:     PetscCall(BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep));

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

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

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

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

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

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

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

349:   /* Get the starting Arnoldi vector */
350:   PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
351:   PetscCall(DSSetDimensions(pep->ds,1,PETSC_DETERMINE,PETSC_DETERMINE));
352:   PetscCall(BVSetActiveColumns(ctx->V,0,1));
353:   PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
354:   PetscCall(BVGetSignature(ctx->V,vomega));
355:   PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));

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

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

386:     /* Solve projected problem */
387:     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
388:     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
389:     PetscCall(DSUpdateExtraRow(pep->ds));
390:     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));

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

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

414:     /* Update S */
415:     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
416:     PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
417:     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));

419:     /* Copy last column of S */
420:     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
421:     PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
422:     PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DETERMINE,PETSC_DETERMINE));
423:     PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
424:     PetscCall(BVSetSignature(ctx->V,vomega));
425:     PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));

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

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

470:   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
471:   PetscFunctionReturn(PETSC_SUCCESS);
472: }

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

481:   PetscFunctionBegin;
482:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP STOAR Options");

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

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

491:     i = 1;
492:     j = k = PETSC_DECIDE;
493:     PetscCall(PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1));
494:     PetscCall(PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2));
495:     PetscCall(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) PetscCall(PEPSTOARSetDimensions(pep,i,j,k));

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

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

506:   PetscOptionsHeadEnd();
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

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

514:   PetscFunctionBegin;
515:   ctx->lock = lock;
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

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

523:    Logically Collective

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

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

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

538:    Level: advanced

540: .seealso: PEPSTOARGetLocking()
541: @*/
542: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
543: {
544:   PetscFunctionBegin;
547:   PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

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

555:   PetscFunctionBegin;
556:   *lock = ctx->lock;
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }

560: /*@
561:    PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.

563:    Not Collective

565:    Input Parameter:
566: .  pep - the eigenproblem solver context

568:    Output Parameter:
569: .  lock - the locking flag

571:    Level: advanced

573: .seealso: PEPSTOARSetLocking()
574: @*/
575: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
576: {
577:   PetscFunctionBegin;
579:   PetscAssertPointer(lock,2);
580:   PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
581:   PetscFunctionReturn(PETSC_SUCCESS);
582: }

584: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
585: {
586:   PetscInt       i,numsh;
587:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
588:   PEP_SR         sr = ctx->sr;

590:   PetscFunctionBegin;
591:   PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
592:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
593:   switch (pep->state) {
594:   case PEP_STATE_INITIAL:
595:     break;
596:   case PEP_STATE_SETUP:
597:     if (n) *n = 2;
598:     if (shifts) {
599:       PetscCall(PetscMalloc1(2,shifts));
600:       (*shifts)[0] = pep->inta;
601:       (*shifts)[1] = pep->intb;
602:     }
603:     if (inertias) {
604:       PetscCall(PetscMalloc1(2,inertias));
605:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
606:       (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
607:     }
608:     break;
609:   case PEP_STATE_SOLVED:
610:   case PEP_STATE_EIGENVECTORS:
611:     numsh = ctx->nshifts;
612:     if (n) *n = numsh;
613:     if (shifts) {
614:       PetscCall(PetscMalloc1(numsh,shifts));
615:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
616:     }
617:     if (inertias) {
618:       PetscCall(PetscMalloc1(numsh,inertias));
619:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
620:     }
621:     break;
622:   }
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: /*@C
627:    PEPSTOARGetInertias - Gets the values of the shifts and their
628:    corresponding inertias in case of doing spectrum slicing for a
629:    computational interval.

631:    Not Collective

633:    Input Parameter:
634: .  pep - the eigenproblem solver context

636:    Output Parameters:
637: +  n        - number of shifts, including the endpoints of the interval
638: .  shifts   - the values of the shifts used internally in the solver
639: -  inertias - the values of the inertia in each shift

641:    Notes:
642:    If called after PEPSolve(), all shifts used internally by the solver are
643:    returned (including both endpoints and any intermediate ones). If called
644:    before PEPSolve() and after PEPSetUp() then only the information of the
645:    endpoints of subintervals is available.

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

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

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

666:    Level: advanced

668: .seealso: PEPSetInterval()
669: @*/
670: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
671: {
672:   PetscFunctionBegin;
674:   PetscAssertPointer(n,2);
675:   PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
676:   PetscFunctionReturn(PETSC_SUCCESS);
677: }

679: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
680: {
681:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

683:   PetscFunctionBegin;
684:   ctx->detect = detect;
685:   pep->state  = PEP_STATE_INITIAL;
686:   PetscFunctionReturn(PETSC_SUCCESS);
687: }

689: /*@
690:    PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
691:    zeros during the factorizations throughout the spectrum slicing computation.

693:    Logically Collective

695:    Input Parameters:
696: +  pep    - the eigenproblem solver context
697: -  detect - check for zeros

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

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

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

710:    Level: advanced

712: .seealso: PEPSetInterval()
713: @*/
714: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
715: {
716:   PetscFunctionBegin;
719:   PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
720:   PetscFunctionReturn(PETSC_SUCCESS);
721: }

723: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
724: {
725:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

727:   PetscFunctionBegin;
728:   *detect = ctx->detect;
729:   PetscFunctionReturn(PETSC_SUCCESS);
730: }

732: /*@
733:    PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
734:    in spectrum slicing.

736:    Not Collective

738:    Input Parameter:
739: .  pep - the eigenproblem solver context

741:    Output Parameter:
742: .  detect - whether zeros detection is enforced during factorizations

744:    Level: advanced

746: .seealso: PEPSTOARSetDetectZeros()
747: @*/
748: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
749: {
750:   PetscFunctionBegin;
752:   PetscAssertPointer(detect,2);
753:   PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
754:   PetscFunctionReturn(PETSC_SUCCESS);
755: }

757: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
758: {
759:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

761:   PetscFunctionBegin;
762:   PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
763:   ctx->alpha = alpha;
764:   ctx->beta  = beta;
765:   PetscFunctionReturn(PETSC_SUCCESS);
766: }

768: /*@
769:    PEPSTOARSetLinearization - Set the coefficients that define
770:    the linearization of a quadratic eigenproblem.

772:    Logically Collective

774:    Input Parameters:
775: +  pep   - polynomial eigenvalue solver
776: .  alpha - first parameter of the linearization
777: -  beta  - second parameter of the linearization

779:    Options Database Key:
780: .  -pep_stoar_linearization <alpha,beta> - Sets the coefficients

782:    Notes:
783:    Cannot pass zero for both alpha and beta. The default values are
784:    alpha=1 and beta=0.

786:    Level: advanced

788: .seealso: PEPSTOARGetLinearization()
789: @*/
790: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
791: {
792:   PetscFunctionBegin;
796:   PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
797:   PetscFunctionReturn(PETSC_SUCCESS);
798: }

800: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
801: {
802:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

804:   PetscFunctionBegin;
805:   if (alpha) *alpha = ctx->alpha;
806:   if (beta)  *beta  = ctx->beta;
807:   PetscFunctionReturn(PETSC_SUCCESS);
808: }

810: /*@
811:    PEPSTOARGetLinearization - Returns the coefficients that define
812:    the linearization of a quadratic eigenproblem.

814:    Not Collective

816:    Input Parameter:
817: .  pep  - polynomial eigenvalue solver

819:    Output Parameters:
820: +  alpha - the first parameter of the linearization
821: -  beta  - the second parameter of the linearization

823:    Level: advanced

825: .seealso: PEPSTOARSetLinearization()
826: @*/
827: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
828: {
829:   PetscFunctionBegin;
831:   PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
832:   PetscFunctionReturn(PETSC_SUCCESS);
833: }

835: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
836: {
837:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

839:   PetscFunctionBegin;
840:   if (nev != PETSC_CURRENT) {
841:     PetscCheck(nev>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
842:     ctx->nev = nev;
843:   }
844:   if (ncv == PETSC_DETERMINE) {
845:     ctx->ncv = PETSC_DETERMINE;
846:   } else if (ncv != PETSC_CURRENT) {
847:     PetscCheck(ncv>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
848:     ctx->ncv = ncv;
849:   }
850:   if (mpd == PETSC_DETERMINE) {
851:     ctx->mpd = PETSC_DETERMINE;
852:   } else if (mpd != PETSC_CURRENT) {
853:     PetscCheck(mpd>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
854:     ctx->mpd = mpd;
855:   }
856:   pep->state = PEP_STATE_INITIAL;
857:   PetscFunctionReturn(PETSC_SUCCESS);
858: }

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

865:    Logically Collective

867:    Input Parameters:
868: +  pep - the eigenproblem solver context
869: .  nev - number of eigenvalues to compute
870: .  ncv - the maximum dimension of the subspace to be used by the subsolve
871: -  mpd - the maximum dimension allowed for the projected problem

873:    Options Database Key:
874: +  -pep_stoar_nev <nev> - Sets the number of eigenvalues
875: .  -pep_stoar_ncv <ncv> - Sets the dimension of the subspace
876: -  -pep_stoar_mpd <mpd> - Sets the maximum projected dimension

878:    Level: advanced

880: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
881: @*/
882: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
883: {
884:   PetscFunctionBegin;
889:   PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
890:   PetscFunctionReturn(PETSC_SUCCESS);
891: }

893: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
894: {
895:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

897:   PetscFunctionBegin;
898:   if (nev) *nev = ctx->nev;
899:   if (ncv) *ncv = ctx->ncv;
900:   if (mpd) *mpd = ctx->mpd;
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

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

908:    Not Collective

910:    Input Parameter:
911: .  pep - the eigenproblem solver context

913:    Output Parameters:
914: +  nev - number of eigenvalues to compute
915: .  ncv - the maximum dimension of the subspace to be used by the subsolve
916: -  mpd - the maximum dimension allowed for the projected problem

918:    Level: advanced

920: .seealso: PEPSTOARSetDimensions()
921: @*/
922: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
923: {
924:   PetscFunctionBegin;
926:   PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

930: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
931: {
932:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

934:   PetscFunctionBegin;
935:   ctx->checket = checket;
936:   pep->state   = PEP_STATE_INITIAL;
937:   PetscFunctionReturn(PETSC_SUCCESS);
938: }

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

944:    Logically Collective

946:    Input Parameters:
947: +  pep     - the eigenproblem solver context
948: -  checket - check eigenvalue type

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

954:    Notes:
955:    This option is relevant only for spectrum slicing computations, but it is
956:    ignored if the problem type is PEP_HYPERBOLIC.

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

964:    Level: advanced

966: .seealso: PEPSetProblemType(), PEPSetInterval()
967: @*/
968: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
969: {
970:   PetscFunctionBegin;
973:   PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
974:   PetscFunctionReturn(PETSC_SUCCESS);
975: }

977: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
978: {
979:   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;

981:   PetscFunctionBegin;
982:   *checket = ctx->checket;
983:   PetscFunctionReturn(PETSC_SUCCESS);
984: }

986: /*@
987:    PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
988:    check in spectrum slicing.

990:    Not Collective

992:    Input Parameter:
993: .  pep - the eigenproblem solver context

995:    Output Parameter:
996: .  checket - whether eigenvalue type must be checked during spectrum slcing

998:    Level: advanced

1000: .seealso: PEPSTOARSetCheckEigenvalueType()
1001: @*/
1002: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1003: {
1004:   PetscFunctionBegin;
1006:   PetscAssertPointer(checket,2);
1007:   PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1008:   PetscFunctionReturn(PETSC_SUCCESS);
1009: }

1011: static PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1012: {
1013:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
1014:   PetscBool      isascii;

1016:   PetscFunctionBegin;
1017:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1018:   if (isascii) {
1019:     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
1020:     PetscCall(PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
1021:     if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscCall(PetscViewerASCIIPrintf(viewer,"  checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled"));
1022:   }
1023:   PetscFunctionReturn(PETSC_SUCCESS);
1024: }

1026: static PetscErrorCode PEPReset_STOAR(PEP pep)
1027: {
1028:   PetscFunctionBegin;
1029:   if (pep->which==PEP_ALL) PetscCall(PEPReset_STOAR_QSlice(pep));
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: static PetscErrorCode PEPDestroy_STOAR(PEP pep)
1034: {
1035:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;

1037:   PetscFunctionBegin;
1038:   PetscCall(BVDestroy(&ctx->V));
1039:   PetscCall(PetscFree(pep->data));
1040:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL));
1041:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL));
1042:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL));
1043:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL));
1044:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL));
1045:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL));
1046:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL));
1047:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL));
1048:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL));
1049:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL));
1050:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL));
1051:   PetscFunctionReturn(PETSC_SUCCESS);
1052: }

1054: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1055: {
1056:   PEP_STOAR      *ctx;

1058:   PetscFunctionBegin;
1059:   PetscCall(PetscNew(&ctx));
1060:   pep->data = (void*)ctx;

1062:   pep->lineariz = PETSC_TRUE;
1063:   ctx->lock     = PETSC_TRUE;
1064:   ctx->nev      = 1;
1065:   ctx->ncv      = PETSC_DETERMINE;
1066:   ctx->mpd      = PETSC_DETERMINE;
1067:   ctx->alpha    = 1.0;
1068:   ctx->beta     = 0.0;
1069:   ctx->checket  = PETSC_TRUE;

1071:   pep->ops->setup          = PEPSetUp_STOAR;
1072:   pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1073:   pep->ops->destroy        = PEPDestroy_STOAR;
1074:   pep->ops->view           = PEPView_STOAR;
1075:   pep->ops->backtransform  = PEPBackTransform_Default;
1076:   pep->ops->computevectors = PEPComputeVectors_Default;
1077:   pep->ops->extractvectors = PEPExtractVectors_TOAR;
1078:   pep->ops->reset          = PEPReset_STOAR;

1080:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR));
1081:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR));
1082:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR));
1083:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR));
1084:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR));
1085:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR));
1086:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR));
1087:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR));
1088:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR));
1089:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR));
1090:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR));
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }