Actual source code: stoar.c
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,(PetscErrorCodeFn*)MatMult_STOAR));
92: PetscCall(MatShellSetOperation(Bs[j],MATOP_DESTROY,(PetscErrorCodeFn*)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 polynomial eigensolver context
527: - lock - `PETSC_TRUE` if the locking variant must be selected
529: Options Database Key:
530: . -pep_stoar_locking - sets the locking flag
532: Note:
533: The default is to lock converged eigenpairs when the method restarts.
534: This behavior 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: [](ch:pep), `PEPSTOAR`, `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 polynomial eigensolver context
568: Output Parameter:
569: . lock - the locking flag
571: Level: advanced
573: .seealso: [](ch:pep), `PEPSTOAR`, `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 polynomial eigensolver 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, that is, when
648: an interval has been given with `PEPSetInterval()` and `STSINVERT` is set.
649: See more details in section [](#sec:qslice).
651: The returned arrays should be freed by the user. Can pass `NULL` in any of
652: the two arrays if not required.
654: Fortran Notes:
655: The calling sequence from Fortran is
656: .vb
657: PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
658: PetscInt n
659: PetscReal shifts(*)
660: PetscInt inertias(*)
661: .ve
662: The arrays should be at least of length `n`. The value of `n` can be determined
663: by an initial call
664: .vb
665: PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
666: .ve
668: Level: advanced
670: .seealso: [](ch:pep), [](#sec:qslice), `PEPSTOAR`, `PEPSetInterval()`
671: @*/
672: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
673: {
674: PetscFunctionBegin;
676: PetscAssertPointer(n,2);
677: PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
682: {
683: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
685: PetscFunctionBegin;
686: ctx->detect = detect;
687: pep->state = PEP_STATE_INITIAL;
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: /*@
692: PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
693: zeros during the factorizations throughout the spectrum slicing computation.
695: Logically Collective
697: Input Parameters:
698: + pep - the polynomial eigensolver context
699: - detect - check for zeros
701: Options Database Key:
702: . -pep_stoar_detect_zeros - toggle the zero detection
704: Notes:
705: This flag makes sense only for spectrum slicing runs, that is, when
706: an interval has been given with `PEPSetInterval()` and `STSINVERT` is set.
707: See more details in section [](#sec:qslice).
709: A zero in the factorization indicates that a shift coincides with an eigenvalue.
711: This flag is turned off by default, and may be necessary in some cases.
712: This feature currently requires an external package for factorizations
713: with support for zero detection, e.g. MUMPS.
715: Level: advanced
717: .seealso: [](ch:pep), [](#sec:qslice), `PEPSTOAR`, `PEPSetInterval()`
718: @*/
719: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
720: {
721: PetscFunctionBegin;
724: PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
729: {
730: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
732: PetscFunctionBegin;
733: *detect = ctx->detect;
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: /*@
738: PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
739: in spectrum slicing.
741: Not Collective
743: Input Parameter:
744: . pep - the polynomial eigensolver context
746: Output Parameter:
747: . detect - whether zeros detection is enforced during factorizations
749: Level: advanced
751: .seealso: [](ch:pep), `PEPSTOAR`, `PEPSTOARSetDetectZeros()`
752: @*/
753: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
754: {
755: PetscFunctionBegin;
757: PetscAssertPointer(detect,2);
758: PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
763: {
764: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
766: PetscFunctionBegin;
767: PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
768: ctx->alpha = alpha;
769: ctx->beta = beta;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@
774: PEPSTOARSetLinearization - Set the coefficients that define the
775: symmetric linearization of a quadratic eigenproblem used in STOAR.
777: Logically Collective
779: Input Parameters:
780: + pep - the polynomial eigensolver context
781: . alpha - first parameter of the linearization
782: - beta - second parameter of the linearization
784: Options Database Key:
785: . -pep_stoar_linearization <alpha,beta> - sets the coefficients
787: Notes:
788: See section [](#sec:linearization) for the general expression of
789: the symmetric linearization.
791: Cannot pass zero for both `alpha` and `beta`. The default values are
792: `alpha`=1 and `beta`=0.
794: Level: advanced
796: .seealso: [](ch:pep), [](#sec:linearization), `PEPSTOAR`, `PEPSTOARGetLinearization()`
797: @*/
798: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
799: {
800: PetscFunctionBegin;
804: PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
809: {
810: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
812: PetscFunctionBegin;
813: if (alpha) *alpha = ctx->alpha;
814: if (beta) *beta = ctx->beta;
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: /*@
819: PEPSTOARGetLinearization - Returns the coefficients that define
820: the linearization of a quadratic eigenproblem.
822: Not Collective
824: Input Parameter:
825: . pep - the polynomial eigensolver context
827: Output Parameters:
828: + alpha - the first parameter of the linearization
829: - beta - the second parameter of the linearization
831: Level: advanced
833: .seealso: [](ch:pep), `PEPSTOAR`, `PEPSTOARSetLinearization()`
834: @*/
835: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
836: {
837: PetscFunctionBegin;
839: PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
844: {
845: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
847: PetscFunctionBegin;
848: if (nev != PETSC_CURRENT) {
849: PetscCheck(nev>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
850: ctx->nev = nev;
851: }
852: if (ncv == PETSC_DETERMINE) {
853: ctx->ncv = PETSC_DETERMINE;
854: } else if (ncv != PETSC_CURRENT) {
855: PetscCheck(ncv>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
856: ctx->ncv = ncv;
857: }
858: if (mpd == PETSC_DETERMINE) {
859: ctx->mpd = PETSC_DETERMINE;
860: } else if (mpd != PETSC_CURRENT) {
861: PetscCheck(mpd>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
862: ctx->mpd = mpd;
863: }
864: pep->state = PEP_STATE_INITIAL;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@
869: PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
870: step in case of doing spectrum slicing for a computational interval.
872: Logically Collective
874: Input Parameters:
875: + pep - the polynomial eigensolver context
876: . nev - number of eigenvalues to compute
877: . ncv - the maximum dimension of the subspace to be used by the subsolve
878: - mpd - the maximum dimension allowed for the projected problem
880: Options Database Keys:
881: + -pep_stoar_nev \<nev\> - sets the number of eigenvalues
882: . -pep_stoar_ncv \<ncv\> - sets the dimension of the subspace
883: - -pep_stoar_mpd \<mpd\> - sets the maximum projected dimension
885: Notes:
886: These parameters are relevant only for spectrum slicing runs, that is, when
887: an interval has been given with `PEPSetInterval()` and `STSINVERT` is set.
888: See more details in section [](#sec:qslice).
890: The meaning of the parameters is the same as in `PEPSetDimensions()`, but
891: the ones here apply to every subsolve done by the child `PEP` object.
893: Use `PETSC_DETERMINE` for `ncv` and `mpd` to assign a default value. For any
894: of the arguments, use `PETSC_CURRENT` to preserve the current value.
896: Level: advanced
898: .seealso: [](ch:pep), [](#sec:qslice), `PEPSTOAR`, `PEPSTOARGetDimensions()`, `PEPSetDimensions()`, `PEPSetInterval()`
899: @*/
900: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
901: {
902: PetscFunctionBegin;
907: PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
912: {
913: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
915: PetscFunctionBegin;
916: if (nev) *nev = ctx->nev;
917: if (ncv) *ncv = ctx->ncv;
918: if (mpd) *mpd = ctx->mpd;
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: /*@
923: PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
924: step in case of doing spectrum slicing for a computational interval.
926: Not Collective
928: Input Parameter:
929: . pep - the polynomial eigensolver context
931: Output Parameters:
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: Level: advanced
938: .seealso: [](ch:pep), `PEPSTOAR`, `PEPSTOARSetDimensions()`
939: @*/
940: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
941: {
942: PetscFunctionBegin;
944: PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
949: {
950: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
952: PetscFunctionBegin;
953: ctx->checket = checket;
954: pep->state = PEP_STATE_INITIAL;
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /*@
959: PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
960: obtained throughout the spectrum slicing computation have the same definite type.
962: Logically Collective
964: Input Parameters:
965: + pep - the polynomial eigensolver context
966: - checket - check eigenvalue type
968: Options Database Key:
969: . -pep_stoar_check_eigenvalue_type - toggles the check of eigenvalue type
971: Notes:
972: This option is relevant only for spectrum slicing computations, but it is
973: ignored if the problem type is `PEP_HYPERBOLIC`.
975: This flag is turned on by default, to guarantee that the computed eigenvalues
976: have the same type (otherwise the computed solution might be wrong). But since
977: the check is computationally quite expensive, the check may be turned off if
978: the user knows for sure that all eigenvalues in the requested interval have
979: the same type.
981: Level: advanced
983: .seealso: [](ch:pep), `PEPSTOAR`, `PEPSetProblemType()`, `PEPSetInterval()`
984: @*/
985: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
986: {
987: PetscFunctionBegin;
990: PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
995: {
996: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
998: PetscFunctionBegin;
999: *checket = ctx->checket;
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
1005: check in spectrum slicing.
1007: Not Collective
1009: Input Parameter:
1010: . pep - the polynomial eigensolver context
1012: Output Parameter:
1013: . checket - whether eigenvalue type must be checked during spectrum slcing
1015: Level: advanced
1017: .seealso: [](ch:pep), `PEPSTOAR`, `PEPSTOARSetCheckEigenvalueType()`
1018: @*/
1019: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1020: {
1021: PetscFunctionBegin;
1023: PetscAssertPointer(checket,2);
1024: PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: static PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1029: {
1030: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1031: PetscBool isascii;
1033: PetscFunctionBegin;
1034: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1035: if (isascii) {
1036: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1037: PetscCall(PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
1038: if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscCall(PetscViewerASCIIPrintf(viewer," checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled"));
1039: }
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: static PetscErrorCode PEPReset_STOAR(PEP pep)
1044: {
1045: PetscFunctionBegin;
1046: if (pep->which==PEP_ALL) PetscCall(PEPReset_STOAR_QSlice(pep));
1047: PetscFunctionReturn(PETSC_SUCCESS);
1048: }
1050: static PetscErrorCode PEPDestroy_STOAR(PEP pep)
1051: {
1052: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1054: PetscFunctionBegin;
1055: PetscCall(BVDestroy(&ctx->V));
1056: PetscCall(PetscFree(pep->data));
1057: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL));
1058: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL));
1059: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL));
1060: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL));
1061: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL));
1062: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL));
1063: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL));
1064: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL));
1065: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL));
1066: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL));
1067: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL));
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }
1071: /*MC
1072: PEPSTOAR - PEPSTOAR = "stoar" - A symmetric variant of TOAR.
1074: Notes:
1075: This solver is available for quadratic eigenproblems only.
1077: It is an alternative to `PEPTOAR` in the case of `PEP_HERMITIAN`
1078: problems. It tries to exploit the symmetry of the matrices by working
1079: with a special linearization where the matrices are Hermitian. However,
1080: this pencil in indefinite, so it uses a pseudo-Lanczos recurrence that
1081: may become numerically unstable. The method is described in {cite:p}`Cam16c`.
1083: The solver incorporates support for interval computation with `PEPSetInterval()`.
1084: Then it will proceed with a spectrum slicing scheme as described in
1085: {cite:p}`Cam20b`. This is particularly useful in the case of `PEP_HYPERBOLIC`
1086: problems.
1088: Level: beginner
1090: .seealso: [](ch:pep), `PEP`, `PEPType`, `PEPSetType()`
1091: M*/
1092: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1093: {
1094: PEP_STOAR *ctx;
1096: PetscFunctionBegin;
1097: PetscCall(PetscNew(&ctx));
1098: pep->data = (void*)ctx;
1100: pep->lineariz = PETSC_TRUE;
1101: ctx->lock = PETSC_TRUE;
1102: ctx->nev = 1;
1103: ctx->ncv = PETSC_DETERMINE;
1104: ctx->mpd = PETSC_DETERMINE;
1105: ctx->alpha = 1.0;
1106: ctx->beta = 0.0;
1107: ctx->checket = PETSC_TRUE;
1109: pep->ops->setup = PEPSetUp_STOAR;
1110: pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1111: pep->ops->destroy = PEPDestroy_STOAR;
1112: pep->ops->view = PEPView_STOAR;
1113: pep->ops->backtransform = PEPBackTransform_Default;
1114: pep->ops->computevectors = PEPComputeVectors_Default;
1115: pep->ops->extractvectors = PEPExtractVectors_TOAR;
1116: pep->ops->reset = PEPReset_STOAR;
1118: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR));
1119: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR));
1120: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR));
1121: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR));
1122: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR));
1123: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR));
1124: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR));
1125: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR));
1126: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR));
1127: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR));
1128: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR));
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }