Actual source code: stoar.c
slepc-3.17.2 2022-08-09
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc polynomial eigensolver: "stoar"
13: Method: S-TOAR
15: Algorithm:
17: Symmetric Two-Level Orthogonal Arnoldi.
19: References:
21: [1] C. Campos and J.E. Roman, "Restarted Q-Arnoldi-type methods
22: exploiting symmetry in quadratic eigenvalue problems", BIT
23: Numer. Math. 56(4):1213-1236, 2016.
24: */
26: #include <slepc/private/pepimpl.h>
27: #include "../src/pep/impls/krylov/pepkrylov.h"
28: #include <slepcblaslapack.h>
30: static PetscBool cited = PETSC_FALSE;
31: static const char citation[] =
32: "@Article{slepc-stoar,\n"
33: " author = \"C. Campos and J. E. Roman\",\n"
34: " title = \"Restarted {Q-Arnoldi-type} methods exploiting symmetry in quadratic eigenvalue problems\",\n"
35: " journal = \"{BIT} Numer. Math.\",\n"
36: " volume = \"56\",\n"
37: " number = \"4\",\n"
38: " pages = \"1213--1236\",\n"
39: " year = \"2016,\"\n"
40: " doi = \"https://doi.org/10.1007/s10543-016-0601-5\"\n"
41: "}\n";
43: typedef struct {
44: PetscReal scal[2];
45: Mat A[2];
46: Vec t;
47: } PEP_STOAR_MATSHELL;
49: static PetscErrorCode MatMult_STOAR(Mat A,Vec x,Vec y)
50: {
51: PEP_STOAR_MATSHELL *ctx;
53: MatShellGetContext(A,&ctx);
54: MatMult(ctx->A[0],x,y);
55: VecScale(y,ctx->scal[0]);
56: if (ctx->scal[1]) {
57: MatMult(ctx->A[1],x,ctx->t);
58: VecAXPY(y,ctx->scal[1],ctx->t);
59: }
60: PetscFunctionReturn(0);
61: }
63: static PetscErrorCode MatDestroy_STOAR(Mat A)
64: {
65: PEP_STOAR_MATSHELL *ctx;
67: MatShellGetContext(A,&ctx);
68: VecDestroy(&ctx->t);
69: PetscFree(ctx);
70: PetscFunctionReturn(0);
71: }
73: PetscErrorCode PEPSTOARSetUpInnerMatrix(PEP pep,Mat *B)
74: {
75: Mat pB[4],Bs[3],D[3];
76: PetscInt i,j,n,m;
77: PEP_STOAR_MATSHELL *ctxMat[3];
78: PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
80: for (i=0;i<3;i++) {
81: STGetMatrixTransformed(pep->st,i,&D[i]); /* D[2] = M */
82: }
83: MatGetLocalSize(D[2],&m,&n);
85: for (j=0;j<3;j++) {
86: PetscNew(ctxMat+j);
87: MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&Bs[j]);
88: MatShellSetOperation(Bs[j],MATOP_MULT,(void(*)(void))MatMult_STOAR);
89: MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_STOAR);
90: }
91: for (i=0;i<4;i++) pB[i] = NULL;
92: if (ctx->alpha) {
93: ctxMat[0]->A[0] = D[0]; ctxMat[0]->scal[0] = ctx->alpha; ctxMat[0]->scal[1] = 0.0;
94: ctxMat[2]->A[0] = D[2]; ctxMat[2]->scal[0] = -ctx->alpha*pep->sfactor*pep->sfactor; ctxMat[2]->scal[1] = 0.0;
95: pB[0] = Bs[0]; pB[3] = Bs[2];
96: }
97: if (ctx->beta) {
98: i = (ctx->alpha)?1:0;
99: ctxMat[0]->scal[1] = 0.0;
100: ctxMat[0]->A[i] = D[1]; ctxMat[0]->scal[i] = -ctx->beta*pep->sfactor;
101: ctxMat[1]->A[0] = D[2]; ctxMat[1]->scal[0] = -ctx->beta*pep->sfactor*pep->sfactor; ctxMat[1]->scal[1] = 0.0;
102: pB[0] = Bs[0]; pB[1] = pB[2] = Bs[1];
103: }
104: BVCreateVec(pep->V,&ctxMat[0]->t);
105: MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pB,B);
106: for (j=0;j<3;j++) MatDestroy(&Bs[j]);
107: PetscFunctionReturn(0);
108: }
110: PetscErrorCode PEPSetUp_STOAR(PEP pep)
111: {
112: PetscBool sinv,flg;
113: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
114: PetscInt ld,i;
115: PetscReal eta;
116: BVOrthogType otype;
117: BVOrthogBlockType obtype;
119: PEPCheckHermitian(pep);
120: PEPCheckQuadratic(pep);
121: PEPCheckShiftSinvert(pep);
122: /* spectrum slicing requires special treatment of default values */
123: if (pep->which==PEP_ALL) {
124: pep->ops->solve = PEPSolve_STOAR_QSlice;
125: pep->ops->extractvectors = NULL;
126: pep->ops->setdefaultst = NULL;
127: PEPSetUp_STOAR_QSlice(pep);
128: } else {
129: PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
131: if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
132: pep->ops->solve = PEPSolve_STOAR;
133: ld = pep->ncv+2;
134: DSSetType(pep->ds,DSGHIEP);
135: DSSetCompact(pep->ds,PETSC_TRUE);
136: DSSetExtraRow(pep->ds,PETSC_TRUE);
137: DSAllocate(pep->ds,ld);
138: PEPBasisCoefficients(pep,pep->pbc);
139: STGetTransform(pep->st,&flg);
140: if (!flg) {
141: PetscFree(pep->solvematcoeffs);
142: PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
143: PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
144: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
145: if (sinv) PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
146: else {
147: for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
148: pep->solvematcoeffs[pep->nmat-1] = 1.0;
149: }
150: }
151: }
152: if (!pep->which) PEPSetWhichEigenpairs_Default(pep);
153: PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_REGION);
155: PEPAllocateSolution(pep,2);
156: PEPSetWorkVecs(pep,4);
157: BVDestroy(&ctx->V);
158: BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
159: BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype);
160: BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
161: PetscFunctionReturn(0);
162: }
164: /*
165: Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
166: */
167: static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
168: {
169: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
170: PetscInt i,j,m=*M,l,lock;
171: PetscInt lds,d,ld,offq,nqt,ldds;
172: Vec v=t_[0],t=t_[1],q=t_[2];
173: PetscReal norm,sym=0.0,fro=0.0,*f;
174: PetscScalar *y,*S,*x,sigma;
175: PetscBLASInt j_,one=1;
176: PetscBool lindep,flg,sinvert=PETSC_FALSE;
177: Mat MS;
179: PetscMalloc1(*M,&y);
180: BVGetSizes(pep->V,NULL,NULL,&ld);
181: BVTensorGetDegree(ctx->V,&d);
182: BVGetActiveColumns(pep->V,&lock,&nqt);
183: lds = d*ld;
184: offq = ld;
185: DSGetLeadingDimension(pep->ds,&ldds);
186: *breakdown = PETSC_FALSE; /* ----- */
187: DSGetDimensions(pep->ds,NULL,&l,NULL,NULL);
188: BVSetActiveColumns(ctx->V,0,m);
189: BVSetActiveColumns(pep->V,0,nqt);
190: STGetTransform(pep->st,&flg);
191: if (!flg) {
192: /* spectral transformation handled by the solver */
193: PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,"");
195: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert);
196: STGetShift(pep->st,&sigma);
197: }
198: for (j=k;j<m;j++) {
199: /* apply operator */
200: BVTensorGetFactors(ctx->V,NULL,&MS);
201: MatDenseGetArray(MS,&S);
202: BVGetColumn(pep->V,nqt,&t);
203: BVMultVec(pep->V,1.0,0.0,v,S+j*lds);
204: if (!sinvert) {
205: STMatMult(pep->st,0,v,q);
206: BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
207: STMatMult(pep->st,1,v,t);
208: VecAXPY(q,pep->sfactor,t);
209: if (ctx->beta && ctx->alpha) {
210: STMatMult(pep->st,2,v,t);
211: VecAXPY(q,-pep->sfactor*pep->sfactor*ctx->beta/ctx->alpha,t);
212: }
213: STMatSolve(pep->st,q,t);
214: VecScale(t,-1.0/(pep->sfactor*pep->sfactor));
215: } else {
216: STMatMult(pep->st,1,v,q);
217: STMatMult(pep->st,2,v,t);
218: VecAXPY(q,sigma*pep->sfactor,t);
219: VecScale(q,pep->sfactor);
220: BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
221: STMatMult(pep->st,2,v,t);
222: VecAXPY(q,pep->sfactor*pep->sfactor,t);
223: STMatSolve(pep->st,q,t);
224: VecScale(t,-1.0);
225: }
226: BVRestoreColumn(pep->V,nqt,&t);
228: /* orthogonalize */
229: if (!sinvert) x = S+offq+(j+1)*lds;
230: else x = S+(j+1)*lds;
231: BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep);
233: if (!lindep) {
234: if (!sinvert) *(S+offq+(j+1)*lds+nqt) = norm;
235: else *(S+(j+1)*lds+nqt) = norm;
236: BVScaleColumn(pep->V,nqt,1.0/norm);
237: nqt++;
238: }
239: if (!sinvert) {
240: for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
241: if (ctx->beta && ctx->alpha) {
242: for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+offq+i) -= *(S+(j+1)*lds+i)*ctx->beta/ctx->alpha;
243: }
244: } else for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
245: BVSetActiveColumns(pep->V,0,nqt);
246: MatDenseRestoreArray(MS,&S);
247: BVTensorRestoreFactors(ctx->V,NULL,&MS);
249: /* level-2 orthogonalization */
250: BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep);
251: a[j] = PetscRealPart(y[j]);
252: omega[j+1] = (norm > 0)?1.0:-1.0;
253: BVScaleColumn(ctx->V,j+1,1.0/norm);
254: b[j] = PetscAbsReal(norm);
256: /* check symmetry */
257: DSGetArrayReal(pep->ds,DS_MAT_T,&f);
258: if (j==k) {
259: for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
260: for (i=0;i<l;i++) y[i] = 0.0;
261: }
262: DSRestoreArrayReal(pep->ds,DS_MAT_T,&f);
263: if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
264: PetscBLASIntCast(j,&j_);
265: sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
266: fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
267: if (j>0) fro = SlepcAbs(fro,b[j-1]);
268: if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
269: *symmlost = PETSC_TRUE;
270: *M=j;
271: break;
272: }
273: }
274: BVSetActiveColumns(pep->V,lock,nqt);
275: BVSetActiveColumns(ctx->V,0,*M);
276: PetscFree(y);
277: PetscFunctionReturn(0);
278: }
280: #if 0
281: static PetscErrorCode PEPSTOARpreKConvergence(PEP pep,PetscInt nv,PetscReal *norm,Vec *w)
282: {
283: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
284: PetscBLASInt n_,one=1;
285: PetscInt lds=2*ctx->ld;
286: PetscReal t1,t2;
287: PetscScalar *S=ctx->S;
289: PetscBLASIntCast(nv+2,&n_);
290: t1 = BLASnrm2_(&n_,S+nv*2*ctx->ld,&one);
291: t2 = BLASnrm2_(&n_,S+(nv*2+1)*ctx->ld,&one);
292: *norm = SlepcAbs(t1,t2);
293: BVSetActiveColumns(pep->V,0,nv+2);
294: BVMultVec(pep->V,1.0,0.0,w[1],S+nv*lds);
295: STMatMult(pep->st,0,w[1],w[2]);
296: VecNorm(w[2],NORM_2,&t1);
297: BVMultVec(pep->V,1.0,0.0,w[1],S+ctx->ld+nv*lds);
298: STMatMult(pep->st,2,w[1],w[2]);
299: VecNorm(w[2],NORM_2,&t2);
300: t2 *= pep->sfactor*pep->sfactor;
301: *norm = PetscMax(*norm,SlepcAbs(t1,t2));
302: PetscFunctionReturn(0);
303: }
304: #endif
306: PetscErrorCode PEPSolve_STOAR(PEP pep)
307: {
308: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
309: PetscInt j,k,l,nv=0,ld,ldds,t,nq=0;
310: PetscInt nconv=0,deg=pep->nmat-1;
311: PetscScalar *om,sigma;
312: PetscReal beta,norm=1.0,*omega,*a,*b;
313: PetscBool breakdown,symmlost=PETSC_FALSE,sinv=PETSC_FALSE,falselock=PETSC_TRUE,flg;
314: Mat MQ,A;
315: Vec vomega;
317: PetscCitationsRegister(citation,&cited);
318: PEPSTOARSetUpInnerMatrix(pep,&A);
319: BVSetMatrix(ctx->V,A,PETSC_TRUE);
320: MatDestroy(&A);
321: if (ctx->lock) {
322: /* undocumented option to use a cheaper locking instead of the true locking */
323: PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL);
324: }
325: BVGetSizes(pep->V,NULL,NULL,&ld);
326: STGetShift(pep->st,&sigma);
327: STGetTransform(pep->st,&flg);
328: if (pep->sfactor!=1.0) {
329: if (!flg) {
330: pep->target /= pep->sfactor;
331: RGPushScale(pep->rg,1.0/pep->sfactor);
332: STScaleShift(pep->st,1.0/pep->sfactor);
333: sigma /= pep->sfactor;
334: } else {
335: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
336: pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
337: RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
338: STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
339: }
340: }
341: if (flg) sigma = 0.0;
343: /* Get the starting Arnoldi vector */
344: BVTensorBuildFirstColumn(ctx->V,pep->nini);
345: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
346: VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
347: BVSetActiveColumns(ctx->V,0,1);
348: BVGetSignature(ctx->V,vomega);
349: VecGetArray(vomega,&om);
350: omega[0] = PetscRealPart(om[0]);
351: VecRestoreArray(vomega,&om);
352: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
353: VecDestroy(&vomega);
355: /* Restart loop */
356: l = 0;
357: DSGetLeadingDimension(pep->ds,&ldds);
358: while (pep->reason == PEP_CONVERGED_ITERATING) {
359: pep->its++;
360: DSGetArrayReal(pep->ds,DS_MAT_T,&a);
361: b = a+ldds;
362: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
364: /* Compute an nv-step Lanczos factorization */
365: nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
366: PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work);
367: beta = b[nv-1];
368: if (symmlost && nv==pep->nconv+l) {
369: pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
370: pep->nconv = nconv;
371: if (falselock || !ctx->lock) {
372: BVSetActiveColumns(ctx->V,0,pep->nconv);
373: BVTensorCompress(ctx->V,0);
374: }
375: break;
376: }
377: DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
378: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
379: DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l);
380: if (l==0) DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
381: else DSSetState(pep->ds,DS_STATE_RAW);
383: /* Solve projected problem */
384: DSSolve(pep->ds,pep->eigr,pep->eigi);
385: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
386: DSUpdateExtraRow(pep->ds);
387: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
389: /* Check convergence */
390: /* PEPSTOARpreKConvergence(pep,nv,&norm,pep->work);*/
391: norm = 1.0;
392: DSGetDimensions(pep->ds,NULL,NULL,NULL,&t);
393: PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k);
394: (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
396: /* Update l */
397: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
398: else {
399: l = PetscMax(1,(PetscInt)((nv-k)/2));
400: l = PetscMin(l,t);
401: DSGetTruncateSize(pep->ds,k,t,&l);
402: if (!breakdown) {
403: /* Prepare the Rayleigh quotient for restart */
404: DSTruncate(pep->ds,k+l,PETSC_FALSE);
405: }
406: }
407: nconv = k;
408: if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
409: if (l) PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
411: /* Update S */
412: DSGetMat(pep->ds,DS_MAT_Q,&MQ);
413: BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
414: MatDestroy(&MQ);
416: /* Copy last column of S */
417: BVCopyColumn(ctx->V,nv,k+l);
418: DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
419: VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
420: VecGetArray(vomega,&om);
421: for (j=0;j<k+l;j++) om[j] = omega[j];
422: VecRestoreArray(vomega,&om);
423: BVSetActiveColumns(ctx->V,0,k+l);
424: BVSetSignature(ctx->V,vomega);
425: VecDestroy(&vomega);
426: DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
428: if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
429: /* stop if breakdown */
430: PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta);
431: pep->reason = PEP_DIVERGED_BREAKDOWN;
432: }
433: if (pep->reason != PEP_CONVERGED_ITERATING) l--;
434: BVGetActiveColumns(pep->V,NULL,&nq);
435: if (k+l+deg<=nq) {
436: BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
437: if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-pep->nconv);
438: else BVTensorCompress(ctx->V,0);
439: }
440: pep->nconv = k;
441: PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
442: }
444: if (pep->nconv>0) {
445: BVSetActiveColumns(ctx->V,0,pep->nconv);
446: BVGetActiveColumns(pep->V,NULL,&nq);
447: BVSetActiveColumns(pep->V,0,nq);
448: if (nq>pep->nconv) {
449: BVTensorCompress(ctx->V,pep->nconv);
450: BVSetActiveColumns(pep->V,0,pep->nconv);
451: }
452: }
453: STGetTransform(pep->st,&flg);
454: if (!flg && pep->ops->backtransform) (*pep->ops->backtransform)(pep);
455: if (pep->sfactor!=1.0) {
456: for (j=0;j<pep->nconv;j++) {
457: pep->eigr[j] *= pep->sfactor;
458: pep->eigi[j] *= pep->sfactor;
459: }
460: }
461: /* restore original values */
462: if (!flg) {
463: pep->target *= pep->sfactor;
464: STScaleShift(pep->st,pep->sfactor);
465: } else {
466: STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
467: pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
468: }
469: if (pep->sfactor!=1.0) RGPopScale(pep->rg);
471: DSTruncate(pep->ds,pep->nconv,PETSC_TRUE);
472: PetscFunctionReturn(0);
473: }
475: PetscErrorCode PEPSetFromOptions_STOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
476: {
477: PetscBool flg,lock,b,f1,f2,f3;
478: PetscInt i,j,k;
479: PetscReal array[2]={0,0};
480: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
482: PetscOptionsHead(PetscOptionsObject,"PEP STOAR Options");
484: PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg);
485: if (flg) PEPSTOARSetLocking(pep,lock);
487: b = ctx->detect;
488: PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg);
489: if (flg) PEPSTOARSetDetectZeros(pep,b);
491: i = 1;
492: j = k = PETSC_DECIDE;
493: PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1);
494: PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2);
495: PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3);
496: if (f1 || f2 || f3) PEPSTOARSetDimensions(pep,i,j,k);
498: k = 2;
499: PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg);
500: if (flg) PEPSTOARSetLinearization(pep,array[0],array[1]);
502: b = ctx->checket;
503: PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg);
504: if (flg) PEPSTOARSetCheckEigenvalueType(pep,b);
506: PetscOptionsTail();
507: PetscFunctionReturn(0);
508: }
510: static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
511: {
512: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
514: ctx->lock = lock;
515: PetscFunctionReturn(0);
516: }
518: /*@
519: PEPSTOARSetLocking - Choose between locking and non-locking variants of
520: the STOAR method.
522: Logically Collective on pep
524: Input Parameters:
525: + pep - the eigenproblem solver context
526: - lock - true if the locking variant must be selected
528: Options Database Key:
529: . -pep_stoar_locking - Sets the locking flag
531: Notes:
532: The default is to lock converged eigenpairs when the method restarts.
533: This behaviour can be changed so that all directions are kept in the
534: working subspace even if already converged to working accuracy (the
535: non-locking variant).
537: Level: advanced
539: .seealso: PEPSTOARGetLocking()
540: @*/
541: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
542: {
545: PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
546: PetscFunctionReturn(0);
547: }
549: static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
550: {
551: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
553: *lock = ctx->lock;
554: PetscFunctionReturn(0);
555: }
557: /*@
558: PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.
560: Not Collective
562: Input Parameter:
563: . pep - the eigenproblem solver context
565: Output Parameter:
566: . lock - the locking flag
568: Level: advanced
570: .seealso: PEPSTOARSetLocking()
571: @*/
572: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
573: {
576: PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
577: PetscFunctionReturn(0);
578: }
580: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
581: {
582: PetscInt i,numsh;
583: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
584: PEP_SR sr = ctx->sr;
588: switch (pep->state) {
589: case PEP_STATE_INITIAL:
590: break;
591: case PEP_STATE_SETUP:
592: if (n) *n = 2;
593: if (shifts) {
594: PetscMalloc1(2,shifts);
595: (*shifts)[0] = pep->inta;
596: (*shifts)[1] = pep->intb;
597: }
598: if (inertias) {
599: PetscMalloc1(2,inertias);
600: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
601: (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
602: }
603: break;
604: case PEP_STATE_SOLVED:
605: case PEP_STATE_EIGENVECTORS:
606: numsh = ctx->nshifts;
607: if (n) *n = numsh;
608: if (shifts) {
609: PetscMalloc1(numsh,shifts);
610: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
611: }
612: if (inertias) {
613: PetscMalloc1(numsh,inertias);
614: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
615: }
616: break;
617: }
618: PetscFunctionReturn(0);
619: }
621: /*@C
622: PEPSTOARGetInertias - Gets the values of the shifts and their
623: corresponding inertias in case of doing spectrum slicing for a
624: computational interval.
626: Not Collective
628: Input Parameter:
629: . pep - the eigenproblem solver context
631: Output Parameters:
632: + n - number of shifts, including the endpoints of the interval
633: . shifts - the values of the shifts used internally in the solver
634: - inertias - the values of the inertia in each shift
636: Notes:
637: If called after PEPSolve(), all shifts used internally by the solver are
638: returned (including both endpoints and any intermediate ones). If called
639: before PEPSolve() and after PEPSetUp() then only the information of the
640: endpoints of subintervals is available.
642: This function is only available for spectrum slicing runs.
644: The returned arrays should be freed by the user. Can pass NULL in any of
645: the two arrays if not required.
647: Fortran Notes:
648: The calling sequence from Fortran is
649: .vb
650: PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
651: integer n
652: double precision shifts(*)
653: integer inertias(*)
654: .ve
655: The arrays should be at least of length n. The value of n can be determined
656: by an initial call
657: .vb
658: PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
659: .ve
661: Level: advanced
663: .seealso: PEPSetInterval()
664: @*/
665: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
666: {
669: PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
670: PetscFunctionReturn(0);
671: }
673: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
674: {
675: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
677: ctx->detect = detect;
678: pep->state = PEP_STATE_INITIAL;
679: PetscFunctionReturn(0);
680: }
682: /*@
683: PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
684: zeros during the factorizations throughout the spectrum slicing computation.
686: Logically Collective on pep
688: Input Parameters:
689: + pep - the eigenproblem solver context
690: - detect - check for zeros
692: Options Database Key:
693: . -pep_stoar_detect_zeros - Check for zeros; this takes an optional
694: bool value (0/1/no/yes/true/false)
696: Notes:
697: A zero in the factorization indicates that a shift coincides with an eigenvalue.
699: This flag is turned off by default, and may be necessary in some cases.
700: This feature currently requires an external package for factorizations
701: with support for zero detection, e.g. MUMPS.
703: Level: advanced
705: .seealso: PEPSetInterval()
706: @*/
707: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
708: {
711: PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
712: PetscFunctionReturn(0);
713: }
715: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
716: {
717: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
719: *detect = ctx->detect;
720: PetscFunctionReturn(0);
721: }
723: /*@
724: PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
725: in spectrum slicing.
727: Not Collective
729: Input Parameter:
730: . pep - the eigenproblem solver context
732: Output Parameter:
733: . detect - whether zeros detection is enforced during factorizations
735: Level: advanced
737: .seealso: PEPSTOARSetDetectZeros()
738: @*/
739: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
740: {
743: PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
744: PetscFunctionReturn(0);
745: }
747: static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
748: {
749: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
752: ctx->alpha = alpha;
753: ctx->beta = beta;
754: PetscFunctionReturn(0);
755: }
757: /*@
758: PEPSTOARSetLinearization - Set the coefficients that define
759: the linearization of a quadratic eigenproblem.
761: Logically Collective on pep
763: Input Parameters:
764: + pep - polynomial eigenvalue solver
765: . alpha - first parameter of the linearization
766: - beta - second parameter of the linearization
768: Options Database Key:
769: . -pep_stoar_linearization <alpha,beta> - Sets the coefficients
771: Notes:
772: Cannot pass zero for both alpha and beta. The default values are
773: alpha=1 and beta=0.
775: Level: advanced
777: .seealso: PEPSTOARGetLinearization()
778: @*/
779: PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
780: {
784: PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
785: PetscFunctionReturn(0);
786: }
788: static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
789: {
790: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
792: if (alpha) *alpha = ctx->alpha;
793: if (beta) *beta = ctx->beta;
794: PetscFunctionReturn(0);
795: }
797: /*@
798: PEPSTOARGetLinearization - Returns the coefficients that define
799: the linearization of a quadratic eigenproblem.
801: Not Collective
803: Input Parameter:
804: . pep - polynomial eigenvalue solver
806: Output Parameters:
807: + alpha - the first parameter of the linearization
808: - beta - the second parameter of the linearization
810: Level: advanced
812: .seealso: PEPSTOARSetLinearization()
813: @*/
814: PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
815: {
817: PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
818: PetscFunctionReturn(0);
819: }
821: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
822: {
823: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
826: ctx->nev = nev;
827: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
828: ctx->ncv = PETSC_DEFAULT;
829: } else {
831: ctx->ncv = ncv;
832: }
833: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
834: ctx->mpd = PETSC_DEFAULT;
835: } else {
837: ctx->mpd = mpd;
838: }
839: pep->state = PEP_STATE_INITIAL;
840: PetscFunctionReturn(0);
841: }
843: /*@
844: PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
845: step in case of doing spectrum slicing for a computational interval.
846: The meaning of the parameters is the same as in PEPSetDimensions().
848: Logically Collective on pep
850: Input Parameters:
851: + pep - the eigenproblem solver context
852: . nev - number of eigenvalues to compute
853: . ncv - the maximum dimension of the subspace to be used by the subsolve
854: - mpd - the maximum dimension allowed for the projected problem
856: Options Database Key:
857: + -eps_stoar_nev <nev> - Sets the number of eigenvalues
858: . -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
859: - -eps_stoar_mpd <mpd> - Sets the maximum projected dimension
861: Level: advanced
863: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
864: @*/
865: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
866: {
871: PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
872: PetscFunctionReturn(0);
873: }
875: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
876: {
877: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
879: if (nev) *nev = ctx->nev;
880: if (ncv) *ncv = ctx->ncv;
881: if (mpd) *mpd = ctx->mpd;
882: PetscFunctionReturn(0);
883: }
885: /*@
886: PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
887: step in case of doing spectrum slicing for a computational interval.
889: Not Collective
891: Input Parameter:
892: . pep - the eigenproblem solver context
894: Output Parameters:
895: + nev - number of eigenvalues to compute
896: . ncv - the maximum dimension of the subspace to be used by the subsolve
897: - mpd - the maximum dimension allowed for the projected problem
899: Level: advanced
901: .seealso: PEPSTOARSetDimensions()
902: @*/
903: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
904: {
906: PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
907: PetscFunctionReturn(0);
908: }
910: static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
911: {
912: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
914: ctx->checket = checket;
915: pep->state = PEP_STATE_INITIAL;
916: PetscFunctionReturn(0);
917: }
919: /*@
920: PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
921: obtained throughout the spectrum slicing computation have the same definite type.
923: Logically Collective on pep
925: Input Parameters:
926: + pep - the eigenproblem solver context
927: - checket - check eigenvalue type
929: Options Database Key:
930: . -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
931: bool value (0/1/no/yes/true/false)
933: Notes:
934: This option is relevant only for spectrum slicing computations, but it is
935: ignored if the problem type is PEP_HYPERBOLIC.
937: This flag is turned on by default, to guarantee that the computed eigenvalues
938: have the same type (otherwise the computed solution might be wrong). But since
939: the check is computationally quite expensive, the check may be turned off if
940: the user knows for sure that all eigenvalues in the requested interval have
941: the same type.
943: Level: advanced
945: .seealso: PEPSetProblemType(), PEPSetInterval()
946: @*/
947: PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
948: {
951: PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
952: PetscFunctionReturn(0);
953: }
955: static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
956: {
957: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
959: *checket = ctx->checket;
960: PetscFunctionReturn(0);
961: }
963: /*@
964: PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
965: check in spectrum slicing.
967: Not Collective
969: Input Parameter:
970: . pep - the eigenproblem solver context
972: Output Parameter:
973: . checket - whether eigenvalue type must be checked during spectrum slcing
975: Level: advanced
977: .seealso: PEPSTOARSetCheckEigenvalueType()
978: @*/
979: PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
980: {
983: PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
984: PetscFunctionReturn(0);
985: }
987: PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
988: {
989: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
990: PetscBool isascii;
992: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
993: if (isascii) {
994: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
995: PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta);
996: if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscViewerASCIIPrintf(viewer," checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled");
997: }
998: PetscFunctionReturn(0);
999: }
1001: PetscErrorCode PEPReset_STOAR(PEP pep)
1002: {
1003: if (pep->which==PEP_ALL) PEPReset_STOAR_QSlice(pep);
1004: PetscFunctionReturn(0);
1005: }
1007: PetscErrorCode PEPDestroy_STOAR(PEP pep)
1008: {
1009: PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1011: BVDestroy(&ctx->V);
1012: PetscFree(pep->data);
1013: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL);
1014: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL);
1015: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL);
1016: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL);
1017: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL);
1018: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL);
1019: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL);
1020: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL);
1021: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL);
1022: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL);
1023: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL);
1024: PetscFunctionReturn(0);
1025: }
1027: SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1028: {
1029: PEP_STOAR *ctx;
1031: PetscNewLog(pep,&ctx);
1032: pep->data = (void*)ctx;
1034: pep->lineariz = PETSC_TRUE;
1035: ctx->lock = PETSC_TRUE;
1036: ctx->nev = 1;
1037: ctx->ncv = PETSC_DEFAULT;
1038: ctx->mpd = PETSC_DEFAULT;
1039: ctx->alpha = 1.0;
1040: ctx->beta = 0.0;
1041: ctx->checket = PETSC_TRUE;
1043: pep->ops->setup = PEPSetUp_STOAR;
1044: pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1045: pep->ops->destroy = PEPDestroy_STOAR;
1046: pep->ops->view = PEPView_STOAR;
1047: pep->ops->backtransform = PEPBackTransform_Default;
1048: pep->ops->computevectors = PEPComputeVectors_Default;
1049: pep->ops->extractvectors = PEPExtractVectors_TOAR;
1050: pep->ops->reset = PEPReset_STOAR;
1052: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR);
1053: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR);
1054: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR);
1055: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR);
1056: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR);
1057: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR);
1058: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR);
1059: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR);
1060: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR);
1061: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR);
1062: PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR);
1063: PetscFunctionReturn(0);
1064: }