Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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"
12 :
13 : Method: S-TOAR
14 :
15 : Algorithm:
16 :
17 : Symmetric Two-Level Orthogonal Arnoldi.
18 :
19 : References:
20 :
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 : */
25 :
26 : #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
27 : #include "../src/pep/impls/krylov/pepkrylov.h"
28 : #include <slepcblaslapack.h>
29 :
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";
42 :
43 : typedef struct {
44 : PetscReal scal[2];
45 : Mat A[2];
46 : Vec t;
47 : } PEP_STOAR_MATSHELL;
48 :
49 9064 : static PetscErrorCode MatMult_STOAR(Mat A,Vec x,Vec y)
50 : {
51 9064 : PEP_STOAR_MATSHELL *ctx;
52 :
53 9064 : PetscFunctionBegin;
54 9064 : PetscCall(MatShellGetContext(A,&ctx));
55 9064 : PetscCall(MatMult(ctx->A[0],x,y));
56 9064 : PetscCall(VecScale(y,ctx->scal[0]));
57 9064 : if (ctx->scal[1]) {
58 28 : PetscCall(MatMult(ctx->A[1],x,ctx->t));
59 28 : PetscCall(VecAXPY(y,ctx->scal[1],ctx->t));
60 : }
61 9064 : PetscFunctionReturn(PETSC_SUCCESS);
62 : }
63 :
64 60 : static PetscErrorCode MatDestroy_STOAR(Mat A)
65 : {
66 60 : PEP_STOAR_MATSHELL *ctx;
67 :
68 60 : PetscFunctionBegin;
69 60 : PetscCall(MatShellGetContext(A,&ctx));
70 60 : PetscCall(VecDestroy(&ctx->t));
71 60 : PetscCall(PetscFree(ctx));
72 60 : PetscFunctionReturn(PETSC_SUCCESS);
73 : }
74 :
75 20 : PetscErrorCode PEPSTOARSetUpInnerMatrix(PEP pep,Mat *B)
76 : {
77 20 : Mat pB[4],Bs[3],D[3];
78 20 : PetscInt i,j,n,m;
79 20 : PEP_STOAR_MATSHELL *ctxMat[3];
80 20 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
81 :
82 20 : PetscFunctionBegin;
83 80 : for (i=0;i<3;i++) {
84 60 : PetscCall(STGetMatrixTransformed(pep->st,i,&D[i])); /* D[2] = M */
85 : }
86 20 : PetscCall(MatGetLocalSize(D[2],&m,&n));
87 :
88 80 : for (j=0;j<3;j++) {
89 60 : PetscCall(PetscNew(ctxMat+j));
90 60 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&Bs[j]));
91 60 : PetscCall(MatShellSetOperation(Bs[j],MATOP_MULT,(void(*)(void))MatMult_STOAR));
92 60 : PetscCall(MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_STOAR));
93 : }
94 100 : for (i=0;i<4;i++) pB[i] = NULL;
95 20 : if (ctx->alpha) {
96 19 : ctxMat[0]->A[0] = D[0]; ctxMat[0]->scal[0] = ctx->alpha; ctxMat[0]->scal[1] = 0.0;
97 19 : ctxMat[2]->A[0] = D[2]; ctxMat[2]->scal[0] = -ctx->alpha*pep->sfactor*pep->sfactor; ctxMat[2]->scal[1] = 0.0;
98 19 : pB[0] = Bs[0]; pB[3] = Bs[2];
99 : }
100 20 : if (ctx->beta) {
101 2 : i = (ctx->alpha)?1:0;
102 2 : ctxMat[0]->scal[1] = 0.0;
103 2 : ctxMat[0]->A[i] = D[1]; ctxMat[0]->scal[i] = -ctx->beta*pep->sfactor;
104 2 : ctxMat[1]->A[0] = D[2]; ctxMat[1]->scal[0] = -ctx->beta*pep->sfactor*pep->sfactor; ctxMat[1]->scal[1] = 0.0;
105 2 : pB[0] = Bs[0]; pB[1] = pB[2] = Bs[1];
106 : }
107 20 : PetscCall(BVCreateVec(pep->V,&ctxMat[0]->t));
108 20 : PetscCall(MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pB,B));
109 80 : for (j=0;j<3;j++) PetscCall(MatDestroy(&Bs[j]));
110 20 : PetscFunctionReturn(PETSC_SUCCESS);
111 : }
112 :
113 21 : static PetscErrorCode PEPSetUp_STOAR(PEP pep)
114 : {
115 21 : PetscBool sinv,flg;
116 21 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
117 21 : PetscInt ld,i;
118 21 : PetscReal eta;
119 21 : BVOrthogType otype;
120 21 : BVOrthogBlockType obtype;
121 :
122 21 : PetscFunctionBegin;
123 21 : PEPCheckHermitian(pep);
124 21 : PEPCheckQuadratic(pep);
125 21 : PEPCheckShiftSinvert(pep);
126 : /* spectrum slicing requires special treatment of default values */
127 21 : if (pep->which==PEP_ALL) {
128 8 : pep->ops->solve = PEPSolve_STOAR_QSlice;
129 8 : pep->ops->extractvectors = NULL;
130 8 : pep->ops->setdefaultst = NULL;
131 8 : PetscCall(PEPSetUp_STOAR_QSlice(pep));
132 : } else {
133 13 : PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
134 13 : PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
135 13 : if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
136 13 : pep->ops->solve = PEPSolve_STOAR;
137 13 : ld = pep->ncv+2;
138 13 : PetscCall(DSSetType(pep->ds,DSGHIEP));
139 13 : PetscCall(DSSetCompact(pep->ds,PETSC_TRUE));
140 13 : PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
141 13 : PetscCall(DSAllocate(pep->ds,ld));
142 13 : PetscCall(PEPBasisCoefficients(pep,pep->pbc));
143 13 : PetscCall(STGetTransform(pep->st,&flg));
144 13 : if (!flg) {
145 11 : PetscCall(PetscFree(pep->solvematcoeffs));
146 11 : PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
147 11 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
148 11 : if (sinv) PetscCall(PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL));
149 : else {
150 15 : for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
151 5 : pep->solvematcoeffs[pep->nmat-1] = 1.0;
152 : }
153 : }
154 : }
155 21 : if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
156 21 : PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_REGION);
157 :
158 21 : PetscCall(PEPAllocateSolution(pep,2));
159 21 : PetscCall(PEPSetWorkVecs(pep,4));
160 21 : PetscCall(BVDestroy(&ctx->V));
161 21 : PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
162 21 : PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
163 21 : PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
164 21 : PetscFunctionReturn(PETSC_SUCCESS);
165 : }
166 :
167 : /*
168 : Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
169 : */
170 112 : static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
171 : {
172 112 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
173 112 : PetscInt i,j,m=*M,l,lock;
174 112 : PetscInt lds,d,ld,offq,nqt,ldds;
175 112 : Vec v=t_[0],t=t_[1],q=t_[2];
176 112 : PetscReal norm,sym=0.0,fro=0.0,*f;
177 112 : PetscScalar *y,*S,*x,sigma;
178 112 : PetscBLASInt j_,one=1;
179 112 : PetscBool lindep,flg,sinvert=PETSC_FALSE;
180 112 : Mat MS;
181 :
182 112 : PetscFunctionBegin;
183 112 : PetscCall(PetscMalloc1(*M,&y));
184 112 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
185 112 : PetscCall(BVTensorGetDegree(ctx->V,&d));
186 112 : PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
187 112 : lds = d*ld;
188 112 : offq = ld;
189 112 : PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
190 112 : *breakdown = PETSC_FALSE; /* ----- */
191 112 : PetscCall(DSGetDimensions(pep->ds,NULL,&l,NULL,NULL));
192 112 : PetscCall(BVSetActiveColumns(ctx->V,0,m));
193 112 : PetscCall(BVSetActiveColumns(pep->V,0,nqt));
194 112 : PetscCall(STGetTransform(pep->st,&flg));
195 112 : if (!flg) {
196 : /* spectral transformation handled by the solver */
197 106 : PetscCall(PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,""));
198 106 : PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"ST type not supported for TOAR without transforming matrices");
199 106 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert));
200 106 : PetscCall(STGetShift(pep->st,&sigma));
201 : }
202 1237 : for (j=k;j<m;j++) {
203 : /* apply operator */
204 1125 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
205 1125 : PetscCall(MatDenseGetArray(MS,&S));
206 1125 : PetscCall(BVGetColumn(pep->V,nqt,&t));
207 1125 : PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+j*lds));
208 1125 : if (!sinvert) {
209 834 : PetscCall(STMatMult(pep->st,0,v,q));
210 834 : PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
211 834 : PetscCall(STMatMult(pep->st,1,v,t));
212 834 : PetscCall(VecAXPY(q,pep->sfactor,t));
213 834 : if (ctx->beta && ctx->alpha) {
214 0 : PetscCall(STMatMult(pep->st,2,v,t));
215 0 : PetscCall(VecAXPY(q,-pep->sfactor*pep->sfactor*ctx->beta/ctx->alpha,t));
216 : }
217 834 : PetscCall(STMatSolve(pep->st,q,t));
218 834 : PetscCall(VecScale(t,-1.0/(pep->sfactor*pep->sfactor)));
219 : } else {
220 291 : PetscCall(STMatMult(pep->st,1,v,q));
221 291 : PetscCall(STMatMult(pep->st,2,v,t));
222 291 : PetscCall(VecAXPY(q,sigma*pep->sfactor,t));
223 291 : PetscCall(VecScale(q,pep->sfactor));
224 291 : PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
225 291 : PetscCall(STMatMult(pep->st,2,v,t));
226 291 : PetscCall(VecAXPY(q,pep->sfactor*pep->sfactor,t));
227 291 : PetscCall(STMatSolve(pep->st,q,t));
228 291 : PetscCall(VecScale(t,-1.0));
229 : }
230 1125 : PetscCall(BVRestoreColumn(pep->V,nqt,&t));
231 :
232 : /* orthogonalize */
233 1125 : if (!sinvert) x = S+offq+(j+1)*lds;
234 291 : else x = S+(j+1)*lds;
235 1125 : PetscCall(BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep));
236 :
237 1125 : if (!lindep) {
238 1125 : if (!sinvert) *(S+offq+(j+1)*lds+nqt) = norm;
239 291 : else *(S+(j+1)*lds+nqt) = norm;
240 1125 : PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
241 1125 : nqt++;
242 : }
243 1125 : if (!sinvert) {
244 13719 : for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
245 834 : if (ctx->beta && ctx->alpha) {
246 0 : for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+offq+i) -= *(S+(j+1)*lds+i)*ctx->beta/ctx->alpha;
247 : }
248 4633 : } else for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
249 1125 : PetscCall(BVSetActiveColumns(pep->V,0,nqt));
250 1125 : PetscCall(MatDenseRestoreArray(MS,&S));
251 1125 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
252 :
253 : /* level-2 orthogonalization */
254 1125 : PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
255 1125 : a[j] = PetscRealPart(y[j]);
256 1125 : omega[j+1] = (norm > 0)?1.0:-1.0;
257 1125 : PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
258 1125 : b[j] = PetscAbsReal(norm);
259 :
260 : /* check symmetry */
261 1125 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&f));
262 1125 : if (j==k) {
263 894 : for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
264 137 : for (i=0;i<l;i++) y[i] = 0.0;
265 : }
266 1125 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&f));
267 1125 : if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
268 1125 : PetscCall(PetscBLASIntCast(j,&j_));
269 1125 : sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
270 1125 : fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
271 1125 : if (j>0) fro = SlepcAbs(fro,b[j-1]);
272 1954 : if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
273 0 : *symmlost = PETSC_TRUE;
274 0 : *M=j;
275 0 : break;
276 : }
277 : }
278 112 : PetscCall(BVSetActiveColumns(pep->V,lock,nqt));
279 112 : PetscCall(BVSetActiveColumns(ctx->V,0,*M));
280 112 : PetscCall(PetscFree(y));
281 112 : PetscFunctionReturn(PETSC_SUCCESS);
282 : }
283 :
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;
292 :
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
310 :
311 13 : PetscErrorCode PEPSolve_STOAR(PEP pep)
312 : {
313 13 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
314 13 : PetscInt j,k,l,nv=0,ld,ldds,t,nq=0;
315 13 : PetscInt nconv=0,deg=pep->nmat-1;
316 13 : PetscScalar sigma;
317 13 : PetscReal beta,norm=1.0,*omega,*a,*b;
318 13 : PetscBool breakdown,symmlost=PETSC_FALSE,sinv=PETSC_FALSE,falselock=PETSC_TRUE,flg;
319 13 : Mat MQ,A,D;
320 13 : Vec vomega;
321 :
322 13 : PetscFunctionBegin;
323 13 : PetscCall(PetscCitationsRegister(citation,&cited));
324 13 : PetscCall(PEPSTOARSetUpInnerMatrix(pep,&A));
325 13 : PetscCall(BVSetMatrix(ctx->V,A,PETSC_TRUE));
326 13 : PetscCall(MatDestroy(&A));
327 13 : if (ctx->lock) {
328 : /* undocumented option to use a cheaper locking instead of the true locking */
329 12 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL));
330 : }
331 13 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
332 13 : PetscCall(STGetShift(pep->st,&sigma));
333 13 : PetscCall(STGetTransform(pep->st,&flg));
334 13 : if (pep->sfactor!=1.0) {
335 5 : if (!flg) {
336 4 : pep->target /= pep->sfactor;
337 4 : PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
338 4 : PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
339 4 : sigma /= pep->sfactor;
340 : } else {
341 1 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
342 1 : pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
343 1 : PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
344 1 : PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
345 : }
346 : }
347 13 : if (flg) sigma = 0.0;
348 :
349 : /* Get the starting Arnoldi vector */
350 13 : PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
351 13 : PetscCall(DSSetDimensions(pep->ds,1,PETSC_DETERMINE,PETSC_DETERMINE));
352 13 : PetscCall(BVSetActiveColumns(ctx->V,0,1));
353 13 : PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
354 13 : PetscCall(BVGetSignature(ctx->V,vomega));
355 13 : PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
356 :
357 : /* Restart loop */
358 13 : l = 0;
359 13 : PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
360 125 : while (pep->reason == PEP_CONVERGED_ITERATING) {
361 112 : pep->its++;
362 112 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
363 112 : b = a+ldds;
364 112 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));
365 :
366 : /* Compute an nv-step Lanczos factorization */
367 112 : nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
368 112 : PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
369 112 : PetscCall(PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work));
370 112 : beta = b[nv-1];
371 112 : if (symmlost && nv==pep->nconv+l) {
372 0 : pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
373 0 : pep->nconv = nconv;
374 0 : if (falselock || !ctx->lock) {
375 0 : PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
376 0 : PetscCall(BVTensorCompress(ctx->V,0));
377 : }
378 : break;
379 : }
380 112 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
381 112 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega));
382 112 : PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
383 112 : if (l==0) PetscCall(DSSetState(pep->ds,DS_STATE_INTERMEDIATE));
384 99 : else PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
385 :
386 : /* Solve projected problem */
387 112 : PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
388 112 : PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
389 112 : PetscCall(DSUpdateExtraRow(pep->ds));
390 112 : PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
391 :
392 : /* Check convergence */
393 : /* PetscCall(PEPSTOARpreKConvergence(pep,nv,&norm,pep->work));*/
394 112 : norm = 1.0;
395 112 : PetscCall(DSGetDimensions(pep->ds,NULL,NULL,NULL,&t));
396 112 : PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k));
397 112 : PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
398 :
399 : /* Update l */
400 112 : if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
401 : else {
402 99 : l = PetscMax(1,(PetscInt)((nv-k)/2));
403 99 : l = PetscMin(l,t);
404 99 : PetscCall(DSGetTruncateSize(pep->ds,k,t,&l));
405 99 : if (!breakdown) {
406 : /* Prepare the Rayleigh quotient for restart */
407 99 : PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
408 : }
409 : }
410 112 : nconv = k;
411 112 : if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
412 112 : if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
413 :
414 : /* Update S */
415 112 : PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
416 112 : PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
417 112 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
418 :
419 : /* Copy last column of S */
420 112 : PetscCall(BVCopyColumn(ctx->V,nv,k+l));
421 112 : PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
422 112 : PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DETERMINE,PETSC_DETERMINE));
423 112 : PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
424 112 : PetscCall(BVSetSignature(ctx->V,vomega));
425 112 : PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
426 :
427 112 : if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
428 : /* stop if breakdown */
429 0 : PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
430 0 : pep->reason = PEP_DIVERGED_BREAKDOWN;
431 : }
432 112 : if (pep->reason != PEP_CONVERGED_ITERATING) l--;
433 112 : PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
434 112 : if (k+l+deg<=nq) {
435 112 : PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
436 112 : if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
437 112 : else PetscCall(BVTensorCompress(ctx->V,0));
438 : }
439 112 : pep->nconv = k;
440 125 : PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
441 : }
442 :
443 13 : if (pep->nconv>0) {
444 13 : PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
445 13 : PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
446 13 : PetscCall(BVSetActiveColumns(pep->V,0,nq));
447 13 : if (nq>pep->nconv) {
448 13 : PetscCall(BVTensorCompress(ctx->V,pep->nconv));
449 13 : PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
450 : }
451 : }
452 13 : PetscCall(STGetTransform(pep->st,&flg));
453 13 : if (!flg) PetscTryTypeMethod(pep,backtransform);
454 13 : if (pep->sfactor!=1.0) {
455 21 : for (j=0;j<pep->nconv;j++) {
456 16 : pep->eigr[j] *= pep->sfactor;
457 16 : pep->eigi[j] *= pep->sfactor;
458 : }
459 : }
460 : /* restore original values */
461 13 : if (!flg) {
462 11 : pep->target *= pep->sfactor;
463 11 : PetscCall(STScaleShift(pep->st,pep->sfactor));
464 : } else {
465 2 : PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
466 2 : pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
467 : }
468 13 : if (pep->sfactor!=1.0) PetscCall(RGPopScale(pep->rg));
469 :
470 13 : PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
471 13 : PetscFunctionReturn(PETSC_SUCCESS);
472 : }
473 :
474 20 : static PetscErrorCode PEPSetFromOptions_STOAR(PEP pep,PetscOptionItems *PetscOptionsObject)
475 : {
476 20 : PetscBool flg,lock,b,f1,f2,f3;
477 20 : PetscInt i,j,k;
478 20 : PetscReal array[2]={0,0};
479 20 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
480 :
481 20 : PetscFunctionBegin;
482 20 : PetscOptionsHeadBegin(PetscOptionsObject,"PEP STOAR Options");
483 :
484 20 : PetscCall(PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg));
485 20 : if (flg) PetscCall(PEPSTOARSetLocking(pep,lock));
486 :
487 20 : b = ctx->detect;
488 20 : PetscCall(PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg));
489 20 : if (flg) PetscCall(PEPSTOARSetDetectZeros(pep,b));
490 :
491 20 : i = 1;
492 20 : j = k = PETSC_DECIDE;
493 20 : PetscCall(PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1));
494 20 : PetscCall(PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2));
495 20 : PetscCall(PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3));
496 20 : if (f1 || f2 || f3) PetscCall(PEPSTOARSetDimensions(pep,i,j,k));
497 :
498 20 : k = 2;
499 20 : PetscCall(PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg));
500 20 : if (flg) PetscCall(PEPSTOARSetLinearization(pep,array[0],array[1]));
501 :
502 20 : b = ctx->checket;
503 20 : PetscCall(PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg));
504 20 : if (flg) PetscCall(PEPSTOARSetCheckEigenvalueType(pep,b));
505 :
506 20 : PetscOptionsHeadEnd();
507 20 : PetscFunctionReturn(PETSC_SUCCESS);
508 : }
509 :
510 2 : static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
511 : {
512 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
513 :
514 2 : PetscFunctionBegin;
515 2 : ctx->lock = lock;
516 2 : PetscFunctionReturn(PETSC_SUCCESS);
517 : }
518 :
519 : /*@
520 : PEPSTOARSetLocking - Choose between locking and non-locking variants of
521 : the STOAR method.
522 :
523 : Logically Collective
524 :
525 : Input Parameters:
526 : + pep - the eigenproblem solver context
527 : - lock - true if the locking variant must be selected
528 :
529 : Options Database Key:
530 : . -pep_stoar_locking - Sets the locking flag
531 :
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).
537 :
538 : Level: advanced
539 :
540 : .seealso: PEPSTOARGetLocking()
541 : @*/
542 2 : PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
543 : {
544 2 : PetscFunctionBegin;
545 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
546 6 : PetscValidLogicalCollectiveBool(pep,lock,2);
547 2 : PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
548 2 : PetscFunctionReturn(PETSC_SUCCESS);
549 : }
550 :
551 2 : static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
552 : {
553 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
554 :
555 2 : PetscFunctionBegin;
556 2 : *lock = ctx->lock;
557 2 : PetscFunctionReturn(PETSC_SUCCESS);
558 : }
559 :
560 : /*@
561 : PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.
562 :
563 : Not Collective
564 :
565 : Input Parameter:
566 : . pep - the eigenproblem solver context
567 :
568 : Output Parameter:
569 : . lock - the locking flag
570 :
571 : Level: advanced
572 :
573 : .seealso: PEPSTOARSetLocking()
574 : @*/
575 2 : PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
576 : {
577 2 : PetscFunctionBegin;
578 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
579 2 : PetscAssertPointer(lock,2);
580 2 : PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
581 2 : PetscFunctionReturn(PETSC_SUCCESS);
582 : }
583 :
584 2 : static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
585 : {
586 2 : PetscInt i,numsh;
587 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
588 2 : PEP_SR sr = ctx->sr;
589 :
590 2 : PetscFunctionBegin;
591 2 : PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
592 2 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
593 2 : switch (pep->state) {
594 : case PEP_STATE_INITIAL:
595 : break;
596 2 : case PEP_STATE_SETUP:
597 2 : if (n) *n = 2;
598 2 : if (shifts) {
599 2 : PetscCall(PetscMalloc1(2,shifts));
600 2 : (*shifts)[0] = pep->inta;
601 2 : (*shifts)[1] = pep->intb;
602 : }
603 2 : if (inertias) {
604 2 : PetscCall(PetscMalloc1(2,inertias));
605 2 : (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
606 2 : (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
607 : }
608 : break;
609 0 : case PEP_STATE_SOLVED:
610 : case PEP_STATE_EIGENVECTORS:
611 0 : numsh = ctx->nshifts;
612 0 : if (n) *n = numsh;
613 0 : if (shifts) {
614 0 : PetscCall(PetscMalloc1(numsh,shifts));
615 0 : for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
616 : }
617 0 : if (inertias) {
618 0 : PetscCall(PetscMalloc1(numsh,inertias));
619 0 : for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
620 : }
621 : break;
622 : }
623 2 : PetscFunctionReturn(PETSC_SUCCESS);
624 : }
625 :
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.
630 :
631 : Not Collective
632 :
633 : Input Parameter:
634 : . pep - the eigenproblem solver context
635 :
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
640 :
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.
646 :
647 : This function is only available for spectrum slicing runs.
648 :
649 : The returned arrays should be freed by the user. Can pass NULL in any of
650 : the two arrays if not required.
651 :
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
665 :
666 : Level: advanced
667 :
668 : .seealso: PEPSetInterval()
669 : @*/
670 2 : PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
671 : {
672 2 : PetscFunctionBegin;
673 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
674 2 : PetscAssertPointer(n,2);
675 2 : PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
676 2 : PetscFunctionReturn(PETSC_SUCCESS);
677 : }
678 :
679 1 : static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
680 : {
681 1 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
682 :
683 1 : PetscFunctionBegin;
684 1 : ctx->detect = detect;
685 1 : pep->state = PEP_STATE_INITIAL;
686 1 : PetscFunctionReturn(PETSC_SUCCESS);
687 : }
688 :
689 : /*@
690 : PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
691 : zeros during the factorizations throughout the spectrum slicing computation.
692 :
693 : Logically Collective
694 :
695 : Input Parameters:
696 : + pep - the eigenproblem solver context
697 : - detect - check for zeros
698 :
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)
702 :
703 : Notes:
704 : A zero in the factorization indicates that a shift coincides with an eigenvalue.
705 :
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.
709 :
710 : Level: advanced
711 :
712 : .seealso: PEPSetInterval()
713 : @*/
714 1 : PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
715 : {
716 1 : PetscFunctionBegin;
717 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
718 3 : PetscValidLogicalCollectiveBool(pep,detect,2);
719 1 : PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
720 1 : PetscFunctionReturn(PETSC_SUCCESS);
721 : }
722 :
723 2 : static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
724 : {
725 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
726 :
727 2 : PetscFunctionBegin;
728 2 : *detect = ctx->detect;
729 2 : PetscFunctionReturn(PETSC_SUCCESS);
730 : }
731 :
732 : /*@
733 : PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
734 : in spectrum slicing.
735 :
736 : Not Collective
737 :
738 : Input Parameter:
739 : . pep - the eigenproblem solver context
740 :
741 : Output Parameter:
742 : . detect - whether zeros detection is enforced during factorizations
743 :
744 : Level: advanced
745 :
746 : .seealso: PEPSTOARSetDetectZeros()
747 : @*/
748 2 : PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
749 : {
750 2 : PetscFunctionBegin;
751 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
752 2 : PetscAssertPointer(detect,2);
753 2 : PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
754 2 : PetscFunctionReturn(PETSC_SUCCESS);
755 : }
756 :
757 4 : static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
758 : {
759 4 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
760 :
761 4 : PetscFunctionBegin;
762 4 : PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
763 4 : ctx->alpha = alpha;
764 4 : ctx->beta = beta;
765 4 : PetscFunctionReturn(PETSC_SUCCESS);
766 : }
767 :
768 : /*@
769 : PEPSTOARSetLinearization - Set the coefficients that define
770 : the linearization of a quadratic eigenproblem.
771 :
772 : Logically Collective
773 :
774 : Input Parameters:
775 : + pep - polynomial eigenvalue solver
776 : . alpha - first parameter of the linearization
777 : - beta - second parameter of the linearization
778 :
779 : Options Database Key:
780 : . -pep_stoar_linearization <alpha,beta> - Sets the coefficients
781 :
782 : Notes:
783 : Cannot pass zero for both alpha and beta. The default values are
784 : alpha=1 and beta=0.
785 :
786 : Level: advanced
787 :
788 : .seealso: PEPSTOARGetLinearization()
789 : @*/
790 4 : PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
791 : {
792 4 : PetscFunctionBegin;
793 4 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
794 12 : PetscValidLogicalCollectiveReal(pep,alpha,2);
795 12 : PetscValidLogicalCollectiveReal(pep,beta,3);
796 4 : PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
797 4 : PetscFunctionReturn(PETSC_SUCCESS);
798 : }
799 :
800 1 : static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
801 : {
802 1 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
803 :
804 1 : PetscFunctionBegin;
805 1 : if (alpha) *alpha = ctx->alpha;
806 1 : if (beta) *beta = ctx->beta;
807 1 : PetscFunctionReturn(PETSC_SUCCESS);
808 : }
809 :
810 : /*@
811 : PEPSTOARGetLinearization - Returns the coefficients that define
812 : the linearization of a quadratic eigenproblem.
813 :
814 : Not Collective
815 :
816 : Input Parameter:
817 : . pep - polynomial eigenvalue solver
818 :
819 : Output Parameters:
820 : + alpha - the first parameter of the linearization
821 : - beta - the second parameter of the linearization
822 :
823 : Level: advanced
824 :
825 : .seealso: PEPSTOARSetLinearization()
826 : @*/
827 1 : PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
828 : {
829 1 : PetscFunctionBegin;
830 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
831 1 : PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
832 1 : PetscFunctionReturn(PETSC_SUCCESS);
833 : }
834 :
835 3 : static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
836 : {
837 3 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
838 :
839 3 : PetscFunctionBegin;
840 3 : if (nev != PETSC_CURRENT) {
841 3 : PetscCheck(nev>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
842 3 : ctx->nev = nev;
843 : }
844 3 : if (ncv == PETSC_DETERMINE) {
845 2 : ctx->ncv = PETSC_DETERMINE;
846 1 : } else if (ncv != PETSC_CURRENT) {
847 1 : PetscCheck(ncv>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
848 1 : ctx->ncv = ncv;
849 : }
850 3 : if (mpd == PETSC_DETERMINE) {
851 2 : ctx->mpd = PETSC_DETERMINE;
852 1 : } else if (mpd != PETSC_CURRENT) {
853 1 : PetscCheck(mpd>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
854 1 : ctx->mpd = mpd;
855 : }
856 3 : pep->state = PEP_STATE_INITIAL;
857 3 : PetscFunctionReturn(PETSC_SUCCESS);
858 : }
859 :
860 : /*@
861 : PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
862 : step in case of doing spectrum slicing for a computational interval.
863 : The meaning of the parameters is the same as in PEPSetDimensions().
864 :
865 : Logically Collective
866 :
867 : Input Parameters:
868 : + pep - the eigenproblem solver context
869 : . nev - number of eigenvalues to compute
870 : . ncv - the maximum dimension of the subspace to be used by the subsolve
871 : - mpd - the maximum dimension allowed for the projected problem
872 :
873 : Options Database Key:
874 : + -pep_stoar_nev <nev> - Sets the number of eigenvalues
875 : . -pep_stoar_ncv <ncv> - Sets the dimension of the subspace
876 : - -pep_stoar_mpd <mpd> - Sets the maximum projected dimension
877 :
878 : Level: advanced
879 :
880 : .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
881 : @*/
882 3 : PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
883 : {
884 3 : PetscFunctionBegin;
885 3 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
886 9 : PetscValidLogicalCollectiveInt(pep,nev,2);
887 9 : PetscValidLogicalCollectiveInt(pep,ncv,3);
888 9 : PetscValidLogicalCollectiveInt(pep,mpd,4);
889 3 : PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
890 3 : PetscFunctionReturn(PETSC_SUCCESS);
891 : }
892 :
893 2 : static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
894 : {
895 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
896 :
897 2 : PetscFunctionBegin;
898 2 : if (nev) *nev = ctx->nev;
899 2 : if (ncv) *ncv = ctx->ncv;
900 2 : if (mpd) *mpd = ctx->mpd;
901 2 : PetscFunctionReturn(PETSC_SUCCESS);
902 : }
903 :
904 : /*@
905 : PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
906 : step in case of doing spectrum slicing for a computational interval.
907 :
908 : Not Collective
909 :
910 : Input Parameter:
911 : . pep - the eigenproblem solver context
912 :
913 : Output Parameters:
914 : + nev - number of eigenvalues to compute
915 : . ncv - the maximum dimension of the subspace to be used by the subsolve
916 : - mpd - the maximum dimension allowed for the projected problem
917 :
918 : Level: advanced
919 :
920 : .seealso: PEPSTOARSetDimensions()
921 : @*/
922 2 : PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
923 : {
924 2 : PetscFunctionBegin;
925 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
926 2 : PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
927 2 : PetscFunctionReturn(PETSC_SUCCESS);
928 : }
929 :
930 1 : static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
931 : {
932 1 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
933 :
934 1 : PetscFunctionBegin;
935 1 : ctx->checket = checket;
936 1 : pep->state = PEP_STATE_INITIAL;
937 1 : PetscFunctionReturn(PETSC_SUCCESS);
938 : }
939 :
940 : /*@
941 : PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
942 : obtained throughout the spectrum slicing computation have the same definite type.
943 :
944 : Logically Collective
945 :
946 : Input Parameters:
947 : + pep - the eigenproblem solver context
948 : - checket - check eigenvalue type
949 :
950 : Options Database Key:
951 : . -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
952 : bool value (0/1/no/yes/true/false)
953 :
954 : Notes:
955 : This option is relevant only for spectrum slicing computations, but it is
956 : ignored if the problem type is PEP_HYPERBOLIC.
957 :
958 : This flag is turned on by default, to guarantee that the computed eigenvalues
959 : have the same type (otherwise the computed solution might be wrong). But since
960 : the check is computationally quite expensive, the check may be turned off if
961 : the user knows for sure that all eigenvalues in the requested interval have
962 : the same type.
963 :
964 : Level: advanced
965 :
966 : .seealso: PEPSetProblemType(), PEPSetInterval()
967 : @*/
968 1 : PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
969 : {
970 1 : PetscFunctionBegin;
971 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
972 3 : PetscValidLogicalCollectiveBool(pep,checket,2);
973 1 : PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
974 1 : PetscFunctionReturn(PETSC_SUCCESS);
975 : }
976 :
977 2 : static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
978 : {
979 2 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
980 :
981 2 : PetscFunctionBegin;
982 2 : *checket = ctx->checket;
983 2 : PetscFunctionReturn(PETSC_SUCCESS);
984 : }
985 :
986 : /*@
987 : PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
988 : check in spectrum slicing.
989 :
990 : Not Collective
991 :
992 : Input Parameter:
993 : . pep - the eigenproblem solver context
994 :
995 : Output Parameter:
996 : . checket - whether eigenvalue type must be checked during spectrum slcing
997 :
998 : Level: advanced
999 :
1000 : .seealso: PEPSTOARSetCheckEigenvalueType()
1001 : @*/
1002 2 : PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
1003 : {
1004 2 : PetscFunctionBegin;
1005 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
1006 2 : PetscAssertPointer(checket,2);
1007 2 : PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
1008 2 : PetscFunctionReturn(PETSC_SUCCESS);
1009 : }
1010 :
1011 0 : static PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
1012 : {
1013 0 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1014 0 : PetscBool isascii;
1015 :
1016 0 : PetscFunctionBegin;
1017 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1018 0 : if (isascii) {
1019 0 : PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1020 0 : PetscCall(PetscViewerASCIIPrintf(viewer," linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
1021 0 : if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscCall(PetscViewerASCIIPrintf(viewer," checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled"));
1022 : }
1023 0 : PetscFunctionReturn(PETSC_SUCCESS);
1024 : }
1025 :
1026 21 : static PetscErrorCode PEPReset_STOAR(PEP pep)
1027 : {
1028 21 : PetscFunctionBegin;
1029 21 : if (pep->which==PEP_ALL) PetscCall(PEPReset_STOAR_QSlice(pep));
1030 21 : PetscFunctionReturn(PETSC_SUCCESS);
1031 : }
1032 :
1033 20 : static PetscErrorCode PEPDestroy_STOAR(PEP pep)
1034 : {
1035 20 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1036 :
1037 20 : PetscFunctionBegin;
1038 20 : PetscCall(BVDestroy(&ctx->V));
1039 20 : PetscCall(PetscFree(pep->data));
1040 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL));
1041 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL));
1042 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL));
1043 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL));
1044 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL));
1045 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL));
1046 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL));
1047 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL));
1048 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL));
1049 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL));
1050 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL));
1051 20 : PetscFunctionReturn(PETSC_SUCCESS);
1052 : }
1053 :
1054 20 : SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
1055 : {
1056 20 : PEP_STOAR *ctx;
1057 :
1058 20 : PetscFunctionBegin;
1059 20 : PetscCall(PetscNew(&ctx));
1060 20 : pep->data = (void*)ctx;
1061 :
1062 20 : pep->lineariz = PETSC_TRUE;
1063 20 : ctx->lock = PETSC_TRUE;
1064 20 : ctx->nev = 1;
1065 20 : ctx->ncv = PETSC_DETERMINE;
1066 20 : ctx->mpd = PETSC_DETERMINE;
1067 20 : ctx->alpha = 1.0;
1068 20 : ctx->beta = 0.0;
1069 20 : ctx->checket = PETSC_TRUE;
1070 :
1071 20 : pep->ops->setup = PEPSetUp_STOAR;
1072 20 : pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
1073 20 : pep->ops->destroy = PEPDestroy_STOAR;
1074 20 : pep->ops->view = PEPView_STOAR;
1075 20 : pep->ops->backtransform = PEPBackTransform_Default;
1076 20 : pep->ops->computevectors = PEPComputeVectors_Default;
1077 20 : pep->ops->extractvectors = PEPExtractVectors_TOAR;
1078 20 : pep->ops->reset = PEPReset_STOAR;
1079 :
1080 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR));
1081 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR));
1082 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR));
1083 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR));
1084 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR));
1085 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR));
1086 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR));
1087 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR));
1088 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR));
1089 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR));
1090 20 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR));
1091 20 : PetscFunctionReturn(PETSC_SUCCESS);
1092 : }
|
