Actual source code: stoar.c
slepc-3.21.0 2024-03-30
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_DEFAULT) 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_DEFAULT,PETSC_DEFAULT));
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_DEFAULT,PETSC_DEFAULT));
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: PetscCheck(nev>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
841: ctx->nev = nev;
842: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
843: ctx->ncv = PETSC_DEFAULT;
844: } else {
845: PetscCheck(ncv>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
846: ctx->ncv = ncv;
847: }
848: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
849: ctx->mpd = PETSC_DEFAULT;
850: } else {
851: PetscCheck(mpd>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
852: ctx->mpd = mpd;
853: }
854: pep->state = PEP_STATE_INITIAL;
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*@
859: PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
860: step in case of doing spectrum slicing for a computational interval.
861: The meaning of the parameters is the same as in PEPSetDimensions().
863: Logically Collective
865: Input Parameters:
866: + pep - the eigenproblem solver context
867: . nev - number of eigenvalues to compute
868: . ncv - the maximum dimension of the subspace to be used by the subsolve
869: - mpd - the maximum dimension allowed for the projected problem
871: Options Database Key:
872: + -eps_stoar_nev <nev> - Sets the number of eigenvalues
873: . -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
874: - -eps_stoar_mpd <mpd> - Sets the maximum projected dimension
876: Level: advanced
878: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
879: @*/
880: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
881: {
882: PetscFunctionBegin;
887: PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
892: {
893: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
895: PetscFunctionBegin;
896: if (nev) *nev = ctx->nev;
897: if (ncv) *ncv = ctx->ncv;
898: if (mpd) *mpd = ctx->mpd;
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: /*@
903: PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
904: step in case of doing spectrum slicing for a computational interval.
906: Not Collective
908: Input Parameter:
909: . pep - the eigenproblem solver context
911: Output Parameters:
912: + nev - number of eigenvalues to compute
913: . ncv - the maximum dimension of the subspace to be used by the subsolve
914: - mpd - the maximum dimension allowed for the projected problem
916: Level: advanced
918: .seealso: PEPSTOARSetDimensions()
919: @*/
920: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
921: {
922: PetscFunctionBegin;
924: PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
925: PetscFunctionReturn(PETSC_SUCCESS);
926: }
928: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
929: {
930: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
932: PetscFunctionBegin;
933: ctx->checket = checket;
934: pep->state = PEP_STATE_INITIAL;
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /*@
939: PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
940: obtained throughout the spectrum slicing computation have the same definite type.
942: Logically Collective
944: Input Parameters:
945: + pep - the eigenproblem solver context
946: - checket - check eigenvalue type
948: Options Database Key:
949: . -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
950: bool value (0/1/no/yes/true/false)
952: Notes:
953: This option is relevant only for spectrum slicing computations, but it is
954: ignored if the problem type is PEP_HYPERBOLIC.
956: This flag is turned on by default, to guarantee that the computed eigenvalues
957: have the same type (otherwise the computed solution might be wrong). But since
958: the check is computationally quite expensive, the check may be turned off if
959: the user knows for sure that all eigenvalues in the requested interval have
960: the same type.
962: Level: advanced
964: .seealso: PEPSetProblemType(), PEPSetInterval()
965: @*/
966: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
967: {
968: PetscFunctionBegin;
971: PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
976: {
977: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
979: PetscFunctionBegin;
980: *checket = ctx->checket;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /*@
985: PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
986: check in spectrum slicing.
988: Not Collective
990: Input Parameter:
991: . pep - the eigenproblem solver context
993: Output Parameter:
994: . checket - whether eigenvalue type must be checked during spectrum slcing
996: Level: advanced
998: .seealso: PEPSTOARSetCheckEigenvalueType()
999: @*/
1000: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1001: {
1002: PetscFunctionBegin;
1004: PetscAssertPointer(checket,2);
1005: PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1006: PetscFunctionReturn(PETSC_SUCCESS);
1007: }
1009: static PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1010: {
1011: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1012: PetscBool isascii;
1014: PetscFunctionBegin;
1015: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1016: if (isascii) {
1017: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1018: PetscCall(PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
1019: if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscCall(PetscViewerASCIIPrintf(viewer," checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled"));
1020: }
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: static PetscErrorCode PEPReset_STOAR(PEP pep)
1025: {
1026: PetscFunctionBegin;
1027: if (pep->which==PEP_ALL) PetscCall(PEPReset_STOAR_QSlice(pep));
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: static PetscErrorCode PEPDestroy_STOAR(PEP pep)
1032: {
1033: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1035: PetscFunctionBegin;
1036: PetscCall(BVDestroy(&ctx->V));
1037: PetscCall(PetscFree(pep->data));
1038: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL));
1039: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL));
1040: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL));
1041: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL));
1042: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL));
1043: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL));
1044: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL));
1045: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL));
1046: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL));
1047: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL));
1048: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL));
1049: PetscFunctionReturn(PETSC_SUCCESS);
1050: }
1052: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1053: {
1054: PEP_STOAR *ctx;
1056: PetscFunctionBegin;
1057: PetscCall(PetscNew(&ctx));
1058: pep->data = (void*)ctx;
1060: pep->lineariz = PETSC_TRUE;
1061: ctx->lock = PETSC_TRUE;
1062: ctx->nev = 1;
1063: ctx->ncv = PETSC_DEFAULT;
1064: ctx->mpd = PETSC_DEFAULT;
1065: ctx->alpha = 1.0;
1066: ctx->beta = 0.0;
1067: ctx->checket = PETSC_TRUE;
1069: pep->ops->setup = PEPSetUp_STOAR;
1070: pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1071: pep->ops->destroy = PEPDestroy_STOAR;
1072: pep->ops->view = PEPView_STOAR;
1073: pep->ops->backtransform = PEPBackTransform_Default;
1074: pep->ops->computevectors = PEPComputeVectors_Default;
1075: pep->ops->extractvectors = PEPExtractVectors_TOAR;
1076: pep->ops->reset = PEPReset_STOAR;
1078: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR));
1079: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR));
1080: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR));
1081: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR));
1082: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR));
1083: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR));
1084: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR));
1085: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR));
1086: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR));
1087: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR));
1088: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR));
1089: PetscFunctionReturn(PETSC_SUCCESS);
1090: }