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 with spectrum slicing for symmetric quadratic eigenproblems
14 :
15 : Algorithm:
16 :
17 : Symmetric Two-Level Orthogonal Arnoldi.
18 :
19 : References:
20 :
21 : [1] C. Campos and J.E. Roman, "Inertia-based spectrum slicing
22 : for symmetric quadratic eigenvalue problems", Numer. Linear
23 : Algebra Appl. 27(4):e2293, 2020.
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-slice-qep,\n"
33 : " author = \"C. Campos and J. E. Roman\",\n"
34 : " title = \"Inertia-based spectrum slicing for symmetric quadratic eigenvalue problems\",\n"
35 : " journal = \"Numer. Linear Algebra Appl.\",\n"
36 : " volume = \"27\",\n"
37 : " number = \"4\",\n"
38 : " pages = \"e2293\",\n"
39 : " year = \"2020,\"\n"
40 : " doi = \"https://doi.org/10.1002/nla.2293\"\n"
41 : "}\n";
42 :
43 : #define SLICE_PTOL PETSC_SQRT_MACHINE_EPSILON
44 :
45 16 : static PetscErrorCode PEPQSliceResetSR(PEP pep)
46 : {
47 16 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
48 16 : PEP_SR sr=ctx->sr;
49 16 : PEP_shift s;
50 16 : PetscInt i;
51 :
52 16 : PetscFunctionBegin;
53 16 : if (sr) {
54 : /* Reviewing list of shifts to free memory */
55 8 : s = sr->s0;
56 8 : if (s) {
57 21 : while (s->neighb[1]) {
58 14 : s = s->neighb[1];
59 21 : PetscCall(PetscFree(s->neighb[0]));
60 : }
61 7 : PetscCall(PetscFree(s));
62 : }
63 8 : PetscCall(PetscFree(sr->S));
64 177 : for (i=0;i<pep->nconv;i++) PetscCall(PetscFree(sr->qinfo[i].q));
65 8 : PetscCall(PetscFree(sr->qinfo));
66 32 : for (i=0;i<3;i++) PetscCall(VecDestroy(&sr->v[i]));
67 8 : PetscCall(EPSDestroy(&sr->eps));
68 8 : PetscCall(PetscFree(sr));
69 : }
70 16 : ctx->sr = NULL;
71 16 : PetscFunctionReturn(PETSC_SUCCESS);
72 : }
73 :
74 8 : PetscErrorCode PEPReset_STOAR_QSlice(PEP pep)
75 : {
76 8 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
77 :
78 8 : PetscFunctionBegin;
79 8 : PetscCall(PEPQSliceResetSR(pep));
80 8 : PetscCall(PetscFree(ctx->inertias));
81 8 : PetscCall(PetscFree(ctx->shifts));
82 8 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 : /*
86 : PEPQSliceAllocateSolution - Allocate memory storage for common variables such
87 : as eigenvalues and eigenvectors.
88 : */
89 7 : static PetscErrorCode PEPQSliceAllocateSolution(PEP pep)
90 : {
91 7 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
92 7 : PetscInt k;
93 7 : BVType type;
94 7 : Vec t;
95 7 : PEP_SR sr = ctx->sr;
96 :
97 7 : PetscFunctionBegin;
98 : /* allocate space for eigenvalues and friends */
99 7 : k = PetscMax(1,sr->numEigs);
100 7 : PetscCall(PetscFree4(sr->eigr,sr->eigi,sr->errest,sr->perm));
101 7 : PetscCall(PetscCalloc4(k,&sr->eigr,k,&sr->eigi,k,&sr->errest,k,&sr->perm));
102 7 : PetscCall(PetscFree(sr->qinfo));
103 7 : PetscCall(PetscCalloc1(k,&sr->qinfo));
104 :
105 : /* allocate sr->V and transfer options from pep->V */
106 7 : PetscCall(BVDestroy(&sr->V));
107 7 : PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),&sr->V));
108 7 : if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
109 7 : if (!((PetscObject)pep->V)->type_name) PetscCall(BVSetType(sr->V,BVMAT));
110 : else {
111 7 : PetscCall(BVGetType(pep->V,&type));
112 7 : PetscCall(BVSetType(sr->V,type));
113 : }
114 7 : PetscCall(STMatCreateVecsEmpty(pep->st,&t,NULL));
115 7 : PetscCall(BVSetSizesFromVec(sr->V,t,k+1));
116 7 : PetscCall(VecDestroy(&t));
117 7 : sr->ld = k;
118 7 : PetscCall(PetscFree(sr->S));
119 7 : PetscCall(PetscMalloc1((k+1)*sr->ld*(pep->nmat-1),&sr->S));
120 7 : PetscFunctionReturn(PETSC_SUCCESS);
121 : }
122 :
123 : /* Convergence test to compute positive Ritz values */
124 31 : static PetscErrorCode ConvergedPositive(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
125 : {
126 31 : PetscFunctionBegin;
127 31 : *errest = (PetscRealPart(eigr)>0.0)?0.0:res;
128 31 : PetscFunctionReturn(PETSC_SUCCESS);
129 : }
130 :
131 33 : static PetscErrorCode PEPQSliceMatGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros)
132 : {
133 33 : KSP ksp,kspr;
134 33 : PC pc;
135 33 : Mat F;
136 33 : PetscBool flg;
137 :
138 33 : PetscFunctionBegin;
139 33 : if (!pep->solvematcoeffs) PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
140 33 : if (shift==PETSC_MAX_REAL) { /* Inertia of matrix A[2] */
141 0 : pep->solvematcoeffs[0] = 0.0; pep->solvematcoeffs[1] = 0.0; pep->solvematcoeffs[2] = 1.0;
142 33 : } else PetscCall(PEPEvaluateBasis(pep,shift,0,pep->solvematcoeffs,NULL));
143 33 : PetscCall(STMatSetUp(pep->st,pep->sfactor,pep->solvematcoeffs));
144 33 : PetscCall(STGetKSP(pep->st,&ksp));
145 33 : PetscCall(KSPGetPC(ksp,&pc));
146 33 : PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCREDUNDANT,&flg));
147 33 : if (flg) {
148 0 : PetscCall(PCRedundantGetKSP(pc,&kspr));
149 0 : PetscCall(KSPGetPC(kspr,&pc));
150 : }
151 33 : PetscCall(PCFactorGetMatrix(pc,&F));
152 33 : PetscCall(MatGetInertia(F,inertia,zeros,NULL));
153 33 : PetscFunctionReturn(PETSC_SUCCESS);
154 : }
155 :
156 31 : static PetscErrorCode PEPQSliceGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros,PetscInt correction)
157 : {
158 31 : KSP ksp;
159 31 : Mat P;
160 31 : PetscReal nzshift=0.0,dot;
161 31 : PetscRandom rand;
162 31 : PetscInt nconv;
163 31 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
164 31 : PEP_SR sr=ctx->sr;
165 :
166 31 : PetscFunctionBegin;
167 31 : if (shift >= PETSC_MAX_REAL) { /* Right-open interval */
168 0 : *inertia = 0;
169 31 : } else if (shift <= PETSC_MIN_REAL) {
170 0 : *inertia = 0;
171 0 : if (zeros) *zeros = 0;
172 : } else {
173 : /* If the shift is zero, perturb it to a very small positive value.
174 : The goal is that the nonzero pattern is the same in all cases and reuse
175 : the symbolic factorizations */
176 31 : nzshift = (shift==0.0)? 10.0/PETSC_MAX_REAL: shift;
177 31 : PetscCall(PEPQSliceMatGetInertia(pep,nzshift,inertia,zeros));
178 31 : PetscCall(STSetShift(pep->st,nzshift));
179 : }
180 31 : if (!correction) {
181 17 : if (shift >= PETSC_MAX_REAL) *inertia = 2*pep->n;
182 17 : else if (shift>PETSC_MIN_REAL) {
183 17 : PetscCall(STGetKSP(pep->st,&ksp));
184 17 : PetscCall(KSPGetOperators(ksp,&P,NULL));
185 17 : if (*inertia!=pep->n && !sr->v[0]) {
186 6 : PetscCall(MatCreateVecs(P,&sr->v[0],NULL));
187 6 : PetscCall(VecDuplicate(sr->v[0],&sr->v[1]));
188 6 : PetscCall(VecDuplicate(sr->v[0],&sr->v[2]));
189 6 : PetscCall(BVGetRandomContext(pep->V,&rand));
190 6 : PetscCall(VecSetRandom(sr->v[0],rand));
191 : }
192 17 : if (*inertia<pep->n && *inertia>0) {
193 9 : if (!sr->eps) {
194 5 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&sr->eps));
195 5 : PetscCall(EPSSetProblemType(sr->eps,EPS_HEP));
196 5 : PetscCall(EPSSetWhichEigenpairs(sr->eps,EPS_LARGEST_REAL));
197 : }
198 9 : PetscCall(EPSSetConvergenceTestFunction(sr->eps,ConvergedPositive,NULL,NULL));
199 9 : PetscCall(EPSSetOperators(sr->eps,P,NULL));
200 9 : PetscCall(EPSSolve(sr->eps));
201 9 : PetscCall(EPSGetConverged(sr->eps,&nconv));
202 9 : PetscCheck(nconv,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)nzshift);
203 9 : PetscCall(EPSGetEigenpair(sr->eps,0,NULL,NULL,sr->v[0],sr->v[1]));
204 : }
205 17 : if (*inertia!=pep->n) {
206 12 : PetscCall(MatMult(pep->A[1],sr->v[0],sr->v[1]));
207 12 : PetscCall(MatMult(pep->A[2],sr->v[0],sr->v[2]));
208 12 : PetscCall(VecAXPY(sr->v[1],2*nzshift,sr->v[2]));
209 12 : PetscCall(VecDotRealPart(sr->v[1],sr->v[0],&dot));
210 12 : if (dot>0.0) *inertia = 2*pep->n-*inertia;
211 : }
212 : }
213 14 : } else if (correction<0) *inertia = 2*pep->n-*inertia;
214 31 : PetscFunctionReturn(PETSC_SUCCESS);
215 : }
216 :
217 : /*
218 : Check eigenvalue type - used only in non-hyperbolic problems.
219 : All computed eigenvalues must have the same definite type (positive or negative).
220 : If ini=TRUE the type is available in omega, otherwise we compute an eigenvalue
221 : closest to shift and determine its type.
222 : */
223 279 : static PetscErrorCode PEPQSliceCheckEigenvalueType(PEP pep,PetscReal shift,PetscReal omega,PetscBool ini)
224 : {
225 279 : PEP pep2;
226 279 : ST st;
227 279 : PetscInt nconv;
228 279 : PetscScalar lambda;
229 279 : PetscReal dot;
230 279 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
231 279 : PEP_SR sr=ctx->sr;
232 :
233 279 : PetscFunctionBegin;
234 279 : if (!ini) {
235 275 : PetscCheck(-(omega/(shift*ctx->alpha+ctx->beta))*sr->type>=0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)shift);
236 : } else {
237 4 : PetscCall(PEPCreate(PetscObjectComm((PetscObject)pep),&pep2));
238 4 : PetscCall(PEPSetOptionsPrefix(pep2,((PetscObject)pep)->prefix));
239 4 : PetscCall(PEPAppendOptionsPrefix(pep2,"pep_eigenvalue_type_"));
240 4 : PetscCall(PEPSetTolerances(pep2,PETSC_CURRENT,pep->max_it/4));
241 4 : PetscCall(PEPSetType(pep2,PEPTOAR));
242 4 : PetscCall(PEPSetOperators(pep2,pep->nmat,pep->A));
243 4 : PetscCall(PEPSetWhichEigenpairs(pep2,PEP_TARGET_MAGNITUDE));
244 4 : PetscCall(PEPGetRG(pep2,&pep2->rg));
245 4 : PetscCall(RGSetType(pep2->rg,RGINTERVAL));
246 : #if defined(PETSC_USE_COMPLEX)
247 4 : PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,-PETSC_SQRT_MACHINE_EPSILON,PETSC_SQRT_MACHINE_EPSILON));
248 : #else
249 : PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,0.0,0.0));
250 : #endif
251 4 : pep2->target = shift;
252 4 : st = pep2->st;
253 4 : pep2->st = pep->st;
254 4 : PetscCall(PEPSolve(pep2));
255 4 : PetscCall(PEPGetConverged(pep2,&nconv));
256 4 : if (nconv) {
257 4 : PetscCall(PEPGetEigenpair(pep2,0,&lambda,NULL,pep2->work[0],NULL));
258 4 : PetscCall(MatMult(pep->A[1],pep2->work[0],pep2->work[1]));
259 4 : PetscCall(MatMult(pep->A[2],pep2->work[0],pep2->work[2]));
260 4 : PetscCall(VecAXPY(pep2->work[1],2.0*lambda*pep->sfactor,pep2->work[2]));
261 4 : PetscCall(VecDotRealPart(pep2->work[1],pep2->work[0],&dot));
262 8 : PetscCall(PetscInfo(pep,"lambda=%g, %s type\n",(double)PetscRealPart(lambda),(dot>0.0)?"positive":"negative"));
263 6 : if (!sr->type) sr->type = (dot>0.0)?1:-1;
264 2 : else PetscCheck(sr->type*dot>=0.0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)PetscRealPart(lambda));
265 : }
266 4 : pep2->st = st;
267 4 : PetscCall(PEPDestroy(&pep2));
268 : }
269 279 : PetscFunctionReturn(PETSC_SUCCESS);
270 : }
271 :
272 3 : static inline PetscErrorCode PEPQSliceDiscriminant(PEP pep,Vec u,Vec w,PetscReal *d,PetscReal *smas,PetscReal *smenos)
273 : {
274 3 : PetscReal ap,bp,cp,dis;
275 :
276 3 : PetscFunctionBegin;
277 3 : PetscCall(MatMult(pep->A[0],u,w));
278 3 : PetscCall(VecDotRealPart(w,u,&cp));
279 3 : PetscCall(MatMult(pep->A[1],u,w));
280 3 : PetscCall(VecDotRealPart(w,u,&bp));
281 3 : PetscCall(MatMult(pep->A[2],u,w));
282 3 : PetscCall(VecDotRealPart(w,u,&ap));
283 3 : dis = bp*bp-4*ap*cp;
284 3 : if (dis>=0.0 && smas) {
285 3 : if (ap>0) *smas = (-bp+PetscSqrtReal(dis))/(2*ap);
286 3 : else if (ap<0) *smas = (-bp-PetscSqrtReal(dis))/(2*ap);
287 : else {
288 0 : if (bp >0) *smas = -cp/bp;
289 0 : else *smas = PETSC_MAX_REAL;
290 : }
291 : }
292 3 : if (dis>=0.0 && smenos) {
293 2 : if (ap>0) *smenos = (-bp-PetscSqrtReal(dis))/(2*ap);
294 2 : else if (ap<0) *smenos = (-bp+PetscSqrtReal(dis))/(2*ap);
295 : else {
296 0 : if (bp<0) *smenos = -cp/bp;
297 0 : else *smenos = PETSC_MAX_REAL;
298 : }
299 : }
300 3 : if (d) *d = dis;
301 3 : PetscFunctionReturn(PETSC_SUCCESS);
302 : }
303 :
304 2 : static inline PetscErrorCode PEPQSliceEvaluateQEP(PEP pep,PetscScalar x,Mat M,MatStructure str)
305 : {
306 2 : PetscFunctionBegin;
307 2 : PetscCall(MatCopy(pep->A[0],M,SAME_NONZERO_PATTERN));
308 2 : PetscCall(MatAXPY(M,x,pep->A[1],str));
309 2 : PetscCall(MatAXPY(M,x*x,pep->A[2],str));
310 2 : PetscFunctionReturn(PETSC_SUCCESS);
311 : }
312 :
313 : /*@
314 : PEPCheckDefiniteQEP - Determines if a symmetric/Hermitian quadratic eigenvalue problem
315 : is definite or not.
316 :
317 : Collective
318 :
319 : Input Parameter:
320 : . pep - eigensolver context
321 :
322 : Output Parameters:
323 : + xi - first computed parameter
324 : . mu - second computed parameter
325 : . definite - flag indicating that the problem is definite
326 : - hyperbolic - flag indicating that the problem is hyperbolic
327 :
328 : Notes:
329 : This function is intended for quadratic eigenvalue problems, Q(lambda)=A*lambda^2+B*lambda+C,
330 : with symmetric (or Hermitian) coefficient matrices A,B,C.
331 :
332 : On output, the flag 'definite' may have the values -1 (meaning that the QEP is not
333 : definite), 1 (if the problem is definite), or 0 if the algorithm was not able to
334 : determine whether the problem is definite or not.
335 :
336 : If definite=1, the output flag 'hyperbolic' informs in a similar way about whether the
337 : problem is hyperbolic or not.
338 :
339 : If definite=1, the computed values xi and mu satisfy Q(xi)<0 and Q(mu)>0, as
340 : obtained via the method proposed in [Niendorf and Voss, LAA 2010]. Furthermore, if
341 : hyperbolic=1 then only xi is computed.
342 :
343 : Level: advanced
344 :
345 : .seealso: PEPSetProblemType()
346 : @*/
347 1 : PetscErrorCode PEPCheckDefiniteQEP(PEP pep,PetscReal *xi,PetscReal *mu,PetscInt *definite,PetscInt *hyperbolic)
348 : {
349 1 : PetscRandom rand;
350 1 : Vec u,w;
351 1 : PetscReal d=0.0,s=0.0,sp,mut=0.0,omg=0.0,omgp;
352 1 : PetscInt k,its=10,hyp=0,check=0,nconv,inertia,n;
353 1 : Mat M=NULL;
354 1 : MatStructure str;
355 1 : EPS eps;
356 1 : PetscBool transform,ptypehyp;
357 :
358 1 : PetscFunctionBegin;
359 1 : PetscCheck(pep->problem_type==PEP_HERMITIAN || pep->problem_type==PEP_HYPERBOLIC,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only available for Hermitian (or hyperbolic) problems");
360 1 : ptypehyp = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;
361 1 : if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
362 1 : PetscCall(PEPSetDefaultST(pep));
363 1 : PetscCall(STSetMatrices(pep->st,pep->nmat,pep->A));
364 1 : PetscCall(MatGetSize(pep->A[0],&n,NULL));
365 1 : PetscCall(STGetTransform(pep->st,&transform));
366 1 : PetscCall(STSetTransform(pep->st,PETSC_FALSE));
367 1 : PetscCall(STSetUp(pep->st));
368 1 : PetscCall(MatCreateVecs(pep->A[0],&u,&w));
369 1 : PetscCall(PEPGetBV(pep,&pep->V));
370 1 : PetscCall(BVGetRandomContext(pep->V,&rand));
371 1 : PetscCall(VecSetRandom(u,rand));
372 1 : PetscCall(VecNormalize(u,NULL));
373 1 : PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,NULL));
374 1 : if (d<0.0) check = -1;
375 1 : if (!check) {
376 1 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&eps));
377 1 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
378 1 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
379 1 : PetscCall(EPSSetTolerances(eps,PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON),PETSC_DECIDE));
380 1 : PetscCall(MatDuplicate(pep->A[0],MAT_DO_NOT_COPY_VALUES,&M));
381 1 : PetscCall(STGetMatStructure(pep->st,&str));
382 : }
383 1 : for (k=0;k<its&&!check;k++) {
384 1 : PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
385 1 : PetscCall(EPSSetOperators(eps,M,NULL));
386 1 : PetscCall(EPSSolve(eps));
387 1 : PetscCall(EPSGetConverged(eps,&nconv));
388 1 : if (!nconv) break;
389 1 : PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
390 1 : sp = s;
391 1 : PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
392 2 : if (d<0.0) {check = -1; break;}
393 1 : if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
394 1 : if (s>sp) {hyp = -1;}
395 1 : mut = 2*s-sp;
396 1 : PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
397 1 : if (inertia == n) {check = 1; break;}
398 : }
399 1 : for (;k<its&&!check;k++) {
400 0 : mut = (s-omg)/2;
401 0 : PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
402 0 : if (inertia == n) {check = 1; break;}
403 0 : if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
404 0 : PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
405 0 : PetscCall(EPSSetOperators(eps,M,NULL));
406 0 : PetscCall(EPSSolve(eps));
407 0 : PetscCall(EPSGetConverged(eps,&nconv));
408 0 : if (!nconv) break;
409 0 : PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
410 0 : omgp = omg;
411 0 : PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
412 0 : if (d<0.0) {check = -1; break;}
413 0 : if (omg<omgp) hyp = -1;
414 : }
415 1 : if (check==1) *xi = mut;
416 1 : PetscCheck(hyp!=-1 || !ptypehyp,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Problem does not satisfy hyperbolic test; consider removing the hyperbolicity flag");
417 1 : if (check==1 && hyp==0) {
418 0 : PetscCall(PEPQSliceMatGetInertia(pep,PETSC_MAX_REAL,&inertia,NULL));
419 0 : if (inertia == 0) hyp = 1;
420 0 : else hyp = -1;
421 : }
422 1 : if (check==1 && hyp!=1) {
423 1 : check = 0;
424 1 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
425 1 : for (;k<its&&!check;k++) {
426 1 : PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
427 1 : PetscCall(EPSSetOperators(eps,M,NULL));
428 1 : PetscCall(EPSSolve(eps));
429 1 : PetscCall(EPSGetConverged(eps,&nconv));
430 1 : if (!nconv) break;
431 1 : PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
432 1 : sp = s;
433 1 : PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
434 2 : if (d<0.0) {check = -1; break;}
435 1 : if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
436 1 : mut = 2*s-sp;
437 1 : PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
438 1 : if (inertia == 0) {check = 1; break;}
439 : }
440 1 : for (;k<its&&!check;k++) {
441 0 : mut = (s-omg)/2;
442 0 : PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
443 0 : if (inertia == 0) {check = 1; break;}
444 0 : if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
445 0 : PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
446 0 : PetscCall(EPSSetOperators(eps,M,NULL));
447 0 : PetscCall(EPSSolve(eps));
448 0 : PetscCall(EPSGetConverged(eps,&nconv));
449 0 : if (!nconv) break;
450 0 : PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
451 0 : PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
452 0 : if (d<0.0) {check = -1; break;}
453 : }
454 : }
455 1 : if (check==1) *mu = mut;
456 1 : *definite = check;
457 1 : *hyperbolic = hyp;
458 1 : if (M) PetscCall(MatDestroy(&M));
459 1 : PetscCall(VecDestroy(&u));
460 1 : PetscCall(VecDestroy(&w));
461 1 : PetscCall(EPSDestroy(&eps));
462 1 : PetscCall(STSetTransform(pep->st,transform));
463 1 : PetscFunctionReturn(PETSC_SUCCESS);
464 : }
465 :
466 : /*
467 : Dummy backtransform operation
468 : */
469 0 : static PetscErrorCode PEPBackTransform_Skip(PEP pep)
470 : {
471 0 : PetscFunctionBegin;
472 0 : PetscFunctionReturn(PETSC_SUCCESS);
473 : }
474 :
475 8 : PetscErrorCode PEPSetUp_STOAR_QSlice(PEP pep)
476 : {
477 8 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
478 8 : PEP_SR sr;
479 8 : PetscInt ld,i,zeros=0;
480 8 : SlepcSC sc;
481 8 : PetscReal r;
482 :
483 8 : PetscFunctionBegin;
484 8 : PEPCheckSinvertCayley(pep);
485 8 : PetscCheck(pep->inta<pep->intb,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues unless you provide a computational interval with PEPSetInterval()");
486 8 : PetscCheck(pep->intb<PETSC_MAX_REAL || pep->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
487 8 : PEPCheckUnsupportedCondition(pep,PEP_FEATURE_STOPPING,PETSC_TRUE," (with spectrum slicing)");
488 8 : if (pep->tol==(PetscReal)PETSC_DETERMINE) {
489 : #if defined(PETSC_USE_REAL_SINGLE)
490 : pep->tol = SLEPC_DEFAULT_TOL;
491 : #else
492 : /* use tighter tolerance */
493 8 : pep->tol = SLEPC_DEFAULT_TOL*1e-2;
494 : #endif
495 : }
496 8 : if (ctx->nev==1) ctx->nev = PetscMin(20,pep->n); /* nev not set, use default value */
497 8 : PetscCheck(pep->n<=10 || ctx->nev>=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"nev cannot be less than 10 in spectrum slicing runs");
498 8 : pep->ops->backtransform = PEPBackTransform_Skip;
499 8 : if (pep->max_it==PETSC_DETERMINE) pep->max_it = 100;
500 :
501 : /* create spectrum slicing context and initialize it */
502 8 : PetscCall(PEPQSliceResetSR(pep));
503 8 : PetscCall(PetscNew(&sr));
504 8 : ctx->sr = sr;
505 8 : sr->itsKs = 0;
506 8 : sr->nleap = 0;
507 8 : sr->sPres = NULL;
508 :
509 8 : if (pep->solvematcoeffs) PetscCall(PetscFree(pep->solvematcoeffs));
510 8 : PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
511 8 : if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
512 8 : PetscCall(STSetTransform(pep->st,PETSC_FALSE));
513 8 : PetscCall(STSetUp(pep->st));
514 :
515 8 : ctx->hyperbolic = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;
516 :
517 : /* check presence of ends and finding direction */
518 8 : if (pep->inta > PETSC_MIN_REAL || pep->intb >= PETSC_MAX_REAL) {
519 8 : sr->int0 = pep->inta;
520 8 : sr->int1 = pep->intb;
521 8 : sr->dir = 1;
522 8 : if (pep->intb >= PETSC_MAX_REAL) { /* Right-open interval */
523 0 : sr->hasEnd = PETSC_FALSE;
524 8 : } else sr->hasEnd = PETSC_TRUE;
525 : } else {
526 0 : sr->int0 = pep->intb;
527 0 : sr->int1 = pep->inta;
528 0 : sr->dir = -1;
529 0 : sr->hasEnd = PetscNot(pep->inta <= PETSC_MIN_REAL);
530 : }
531 :
532 : /* compute inertia0 */
533 15 : PetscCall(PEPQSliceGetInertia(pep,sr->int0,&sr->inertia0,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
534 8 : PetscCheck(!zeros || (sr->int0!=pep->inta && sr->int0!=pep->intb),((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in the interval endpoint");
535 8 : if (!ctx->hyperbolic && ctx->checket) PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int0,0.0,PETSC_TRUE));
536 :
537 : /* compute inertia1 */
538 15 : PetscCall(PEPQSliceGetInertia(pep,sr->int1,&sr->inertia1,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
539 8 : PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
540 8 : if (!ctx->hyperbolic && ctx->checket && sr->hasEnd) {
541 2 : PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int1,0.0,PETSC_TRUE));
542 2 : PetscCheck(sr->type || sr->inertia1==sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"No information of eigenvalue type in interval");
543 2 : PetscCheck(!sr->type || sr->inertia1!=sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected");
544 2 : if (sr->dir*(sr->inertia1-sr->inertia0)<0) {
545 0 : sr->intcorr = -1;
546 0 : sr->inertia0 = 2*pep->n-sr->inertia0;
547 0 : sr->inertia1 = 2*pep->n-sr->inertia1;
548 2 : } else sr->intcorr = 1;
549 : } else {
550 6 : if (sr->inertia0<=pep->n && sr->inertia1<=pep->n) sr->intcorr = 1;
551 4 : else if (sr->inertia0>=pep->n && sr->inertia1>=pep->n) sr->intcorr = -1;
552 : }
553 :
554 8 : if (sr->hasEnd) {
555 8 : sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
556 8 : i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
557 : }
558 :
559 : /* number of eigenvalues in interval */
560 8 : sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
561 8 : PetscCall(PetscInfo(pep,"QSlice setup: allocating for %" PetscInt_FMT " eigenvalues in [%g,%g]\n",sr->numEigs,(double)pep->inta,(double)pep->intb));
562 8 : if (sr->numEigs) {
563 7 : PetscCall(PEPQSliceAllocateSolution(pep));
564 7 : PetscCall(PEPSetDimensions_Default(pep,ctx->nev,&ctx->ncv,&ctx->mpd));
565 7 : pep->nev = ctx->nev; pep->ncv = ctx->ncv; pep->mpd = ctx->mpd;
566 7 : ld = ctx->ncv+2;
567 7 : PetscCall(DSSetType(pep->ds,DSGHIEP));
568 7 : PetscCall(DSSetCompact(pep->ds,PETSC_TRUE));
569 7 : PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
570 7 : PetscCall(DSAllocate(pep->ds,ld));
571 7 : PetscCall(DSGetSlepcSC(pep->ds,&sc));
572 7 : sc->rg = NULL;
573 7 : sc->comparison = SlepcCompareLargestMagnitude;
574 7 : sc->comparisonctx = NULL;
575 7 : sc->map = NULL;
576 7 : sc->mapobj = NULL;
577 1 : } else {pep->ncv = 0; pep->nev = 0; pep->mpd = 0;}
578 8 : PetscFunctionReturn(PETSC_SUCCESS);
579 : }
580 :
581 : /*
582 : Fills the fields of a shift structure
583 : */
584 23 : static PetscErrorCode PEPCreateShift(PEP pep,PetscReal val,PEP_shift neighb0,PEP_shift neighb1)
585 : {
586 23 : PEP_shift s,*pending2;
587 23 : PetscInt i;
588 23 : PEP_SR sr;
589 23 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
590 :
591 23 : PetscFunctionBegin;
592 23 : sr = ctx->sr;
593 23 : PetscCall(PetscNew(&s));
594 23 : s->value = val;
595 23 : s->neighb[0] = neighb0;
596 23 : if (neighb0) neighb0->neighb[1] = s;
597 23 : s->neighb[1] = neighb1;
598 23 : if (neighb1) neighb1->neighb[0] = s;
599 23 : s->comp[0] = PETSC_FALSE;
600 23 : s->comp[1] = PETSC_FALSE;
601 23 : s->index = -1;
602 23 : s->neigs = 0;
603 23 : s->nconv[0] = s->nconv[1] = 0;
604 23 : s->nsch[0] = s->nsch[1]=0;
605 : /* Inserts in the stack of pending shifts */
606 : /* If needed, the array is resized */
607 23 : if (sr->nPend >= sr->maxPend) {
608 0 : sr->maxPend *= 2;
609 0 : PetscCall(PetscMalloc1(sr->maxPend,&pending2));
610 0 : for (i=0;i<sr->nPend;i++) pending2[i] = sr->pending[i];
611 0 : PetscCall(PetscFree(sr->pending));
612 0 : sr->pending = pending2;
613 : }
614 23 : sr->pending[sr->nPend++]=s;
615 23 : PetscFunctionReturn(PETSC_SUCCESS);
616 : }
617 :
618 : /* Provides next shift to be computed */
619 22 : static PetscErrorCode PEPExtractShift(PEP pep)
620 : {
621 22 : PetscInt iner,zeros=0;
622 22 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
623 22 : PEP_SR sr;
624 22 : PetscReal newShift,aux;
625 22 : PEP_shift sPres;
626 :
627 22 : PetscFunctionBegin;
628 22 : sr = ctx->sr;
629 22 : if (sr->nPend > 0) {
630 15 : if (sr->dirch) {
631 1 : aux = sr->int1; sr->int1 = sr->int0; sr->int0 = aux;
632 1 : iner = sr->inertia1; sr->inertia1 = sr->inertia0; sr->inertia0 = iner;
633 1 : sr->dir *= -1;
634 1 : PetscCall(PetscFree(sr->s0->neighb[1]));
635 1 : PetscCall(PetscFree(sr->s0));
636 1 : sr->nPend--;
637 1 : PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
638 1 : sr->sPrev = NULL;
639 1 : sr->sPres = sr->pending[--sr->nPend];
640 1 : pep->target = sr->sPres->value;
641 1 : sr->s0 = sr->sPres;
642 1 : pep->reason = PEP_CONVERGED_ITERATING;
643 : } else {
644 14 : sr->sPrev = sr->sPres;
645 14 : sr->sPres = sr->pending[--sr->nPend];
646 : }
647 15 : sPres = sr->sPres;
648 30 : PetscCall(PEPQSliceGetInertia(pep,sPres->value,&iner,ctx->detect?&zeros:NULL,sr->intcorr));
649 15 : if (zeros) {
650 0 : newShift = sPres->value*(1.0+SLICE_PTOL);
651 0 : if (sr->dir*(sPres->neighb[0] && newShift-sPres->neighb[0]->value) < 0) newShift = (sPres->value+sPres->neighb[0]->value)/2;
652 0 : else if (sPres->neighb[1] && sr->dir*(sPres->neighb[1]->value-newShift) < 0) newShift = (sPres->value+sPres->neighb[1]->value)/2;
653 0 : PetscCall(PEPQSliceGetInertia(pep,newShift,&iner,&zeros,sr->intcorr));
654 0 : PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)newShift);
655 0 : sPres->value = newShift;
656 : }
657 15 : sr->sPres->inertia = iner;
658 15 : pep->target = sr->sPres->value;
659 15 : pep->reason = PEP_CONVERGED_ITERATING;
660 15 : pep->its = 0;
661 7 : } else sr->sPres = NULL;
662 22 : PetscFunctionReturn(PETSC_SUCCESS);
663 : }
664 :
665 : /*
666 : Obtains value of subsequent shift
667 : */
668 15 : static PetscErrorCode PEPGetNewShiftValue(PEP pep,PetscInt side,PetscReal *newS)
669 : {
670 15 : PetscReal lambda,d_prev;
671 15 : PetscInt i,idxP;
672 15 : PEP_SR sr;
673 15 : PEP_shift sPres,s;
674 15 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
675 :
676 15 : PetscFunctionBegin;
677 15 : sr = ctx->sr;
678 15 : sPres = sr->sPres;
679 15 : if (sPres->neighb[side]) {
680 : /* Completing a previous interval */
681 3 : if (!sPres->neighb[side]->neighb[side] && sPres->neighb[side]->nconv[side]==0) { /* One of the ends might be too far from eigenvalues */
682 0 : if (side) *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[sr->indexEig-1]]))/2;
683 0 : else *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[0]]))/2;
684 3 : } else *newS=(sPres->value + sPres->neighb[side]->value)/2;
685 : } else { /* (Only for side=1). Creating a new interval. */
686 12 : if (sPres->neigs==0) {/* No value has been accepted*/
687 5 : if (sPres->neighb[0]) {
688 : /* Multiplying by 10 the previous distance */
689 4 : *newS = sPres->value + 10*sr->dir*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
690 4 : sr->nleap++;
691 : /* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
692 4 : PetscCheck(sr->hasEnd || sr->nleap<=5,PetscObjectComm((PetscObject)pep),PETSC_ERR_CONV_FAILED,"Unable to compute the wanted eigenvalues with open interval");
693 : } else { /* First shift */
694 1 : if (pep->nconv != 0) {
695 : /* Unaccepted values give information for next shift */
696 : idxP=0;/* Number of values left from shift */
697 0 : for (i=0;i<pep->nconv;i++) {
698 0 : lambda = PetscRealPart(pep->eigr[i]);
699 0 : if (sr->dir*(lambda - sPres->value) <0) idxP++;
700 : else break;
701 : }
702 : /* Avoiding subtraction of eigenvalues (might be the same).*/
703 0 : if (idxP>0) {
704 0 : d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[0]))/(idxP+0.3);
705 : } else {
706 0 : d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[pep->nconv-1]))/(pep->nconv+0.3);
707 : }
708 0 : *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
709 0 : sr->dirch = PETSC_FALSE;
710 : } else { /* No values found, no information for next shift */
711 1 : PetscCheck(!sr->dirch,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"First shift renders no information");
712 1 : sr->dirch = PETSC_TRUE;
713 1 : *newS = sr->int1;
714 : }
715 : }
716 : } else { /* Accepted values found */
717 7 : sr->dirch = PETSC_FALSE;
718 7 : sr->nleap = 0;
719 : /* Average distance of values in previous subinterval */
720 7 : s = sPres->neighb[0];
721 7 : while (s && PetscAbs(s->inertia - sPres->inertia)==0) {
722 0 : s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
723 : }
724 7 : if (s) {
725 3 : d_prev = PetscAbsReal((sPres->value - s->value)/(sPres->inertia - s->inertia));
726 : } else { /* First shift. Average distance obtained with values in this shift */
727 : /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
728 4 : if (sr->dir*(PetscRealPart(sr->eigr[0])-sPres->value)>0 && PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0]))/PetscRealPart(sr->eigr[0])) > PetscSqrtReal(pep->tol)) {
729 3 : d_prev = PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0])))/(sPres->neigs+0.3);
730 : } else {
731 1 : d_prev = PetscAbsReal(PetscRealPart(sr->eigr[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
732 : }
733 : }
734 : /* Average distance is used for next shift by adding it to value on the right or to shift */
735 7 : if (sr->dir*(PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1]) - sPres->value)>0) {
736 7 : *newS = PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1])+ (sr->dir*d_prev*pep->nev)/2;
737 : } else { /* Last accepted value is on the left of shift. Adding to shift */
738 0 : *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
739 : }
740 : }
741 : /* End of interval can not be surpassed */
742 12 : if (sr->dir*(sr->int1 - *newS) < 0) *newS = sr->int1;
743 : }/* of neighb[side]==null */
744 15 : PetscFunctionReturn(PETSC_SUCCESS);
745 : }
746 :
747 : /*
748 : Function for sorting an array of real values
749 : */
750 43 : static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
751 : {
752 43 : PetscReal re;
753 43 : PetscInt i,j,tmp;
754 :
755 43 : PetscFunctionBegin;
756 1006 : if (!prev) for (i=0;i<nr;i++) perm[i] = i;
757 : /* Insertion sort */
758 964 : for (i=1;i<nr;i++) {
759 921 : re = PetscRealPart(r[perm[i]]);
760 921 : j = i-1;
761 3496 : while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
762 2575 : tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
763 : }
764 : }
765 43 : PetscFunctionReturn(PETSC_SUCCESS);
766 : }
767 :
768 : /* Stores the pairs obtained since the last shift in the global arrays */
769 22 : static PetscErrorCode PEPStoreEigenpairs(PEP pep)
770 : {
771 22 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
772 22 : PetscReal lambda,err,*errest;
773 22 : PetscInt i,*aux,count=0,ndef,ld,nconv=pep->nconv,d=pep->nmat-1,idx;
774 22 : PetscBool iscayley,divide=PETSC_FALSE;
775 22 : PEP_SR sr = ctx->sr;
776 22 : PEP_shift sPres;
777 22 : Vec w,vomega;
778 22 : Mat MS;
779 22 : BV tV;
780 22 : PetscScalar *S,*eigr,*tS,*omega;
781 :
782 22 : PetscFunctionBegin;
783 22 : sPres = sr->sPres;
784 22 : sPres->index = sr->indexEig;
785 :
786 22 : if (nconv>sr->ndef0+sr->ndef1) {
787 : /* Back-transform */
788 21 : PetscCall(STBackTransform(pep->st,nconv,pep->eigr,pep->eigi));
789 450 : for (i=0;i<nconv;i++) {
790 : #if defined(PETSC_USE_COMPLEX)
791 429 : if (PetscImaginaryPart(pep->eigr[i])) pep->eigr[i] = sr->int0-sr->dir;
792 : #else
793 : if (pep->eigi[i]) pep->eigr[i] = sr->int0-sr->dir;
794 : #endif
795 : }
796 21 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STCAYLEY,&iscayley));
797 : /* Sort eigenvalues */
798 21 : PetscCall(sortRealEigenvalues(pep->eigr,pep->perm,nconv,PETSC_FALSE,sr->dir));
799 21 : PetscCall(VecCreateSeq(PETSC_COMM_SELF,nconv,&vomega));
800 21 : PetscCall(BVGetSignature(ctx->V,vomega));
801 21 : PetscCall(VecGetArray(vomega,&omega));
802 21 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
803 21 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
804 21 : PetscCall(MatDenseGetArray(MS,&S));
805 : /* Values stored in global array */
806 21 : PetscCall(PetscCalloc4(nconv,&eigr,nconv,&errest,nconv*nconv*d,&tS,nconv,&aux));
807 21 : ndef = sr->ndef0+sr->ndef1;
808 450 : for (i=0;i<nconv;i++) {
809 429 : lambda = PetscRealPart(pep->eigr[pep->perm[i]]);
810 429 : err = pep->errest[pep->perm[i]];
811 429 : if (sr->dir*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0) {/* Valid value */
812 275 : PetscCheck(sr->indexEig+count-ndef<sr->numEigs,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Unexpected error in Spectrum Slicing");
813 275 : PetscCall(PEPQSliceCheckEigenvalueType(pep,lambda,PetscRealPart(omega[pep->perm[i]]),PETSC_FALSE));
814 275 : eigr[count] = lambda;
815 275 : errest[count] = err;
816 275 : if ((sr->dir*(sPres->value - lambda) > 0) && (sr->dir*(lambda - sPres->ext[0]) > 0)) sPres->nconv[0]++;
817 275 : if ((sr->dir*(lambda - sPres->value) > 0) && (sr->dir*(sPres->ext[1] - lambda) > 0)) sPres->nconv[1]++;
818 275 : PetscCall(PetscArraycpy(tS+count*(d*nconv),S+pep->perm[i]*(d*ld),nconv));
819 275 : PetscCall(PetscArraycpy(tS+count*(d*nconv)+nconv,S+pep->perm[i]*(d*ld)+ld,nconv));
820 275 : count++;
821 : }
822 : }
823 21 : PetscCall(VecRestoreArray(vomega,&omega));
824 21 : PetscCall(VecDestroy(&vomega));
825 296 : for (i=0;i<count;i++) {
826 275 : PetscCall(PetscArraycpy(S+i*(d*ld),tS+i*nconv*d,nconv));
827 275 : PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+i*nconv*d+nconv,nconv));
828 : }
829 21 : PetscCall(MatDenseRestoreArray(MS,&S));
830 21 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
831 21 : PetscCall(BVSetActiveColumns(ctx->V,0,count));
832 21 : PetscCall(BVTensorCompress(ctx->V,count));
833 21 : if (sr->sPres->nconv[0] && sr->sPres->nconv[1]) {
834 6 : divide = PETSC_TRUE;
835 6 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
836 6 : PetscCall(MatDenseGetArray(MS,&S));
837 6 : PetscCall(PetscArrayzero(tS,nconv*nconv*d));
838 137 : for (i=0;i<count;i++) {
839 131 : PetscCall(PetscArraycpy(tS+i*nconv*d,S+i*(d*ld),count));
840 131 : PetscCall(PetscArraycpy(tS+i*nconv*d+nconv,S+i*(d*ld)+ld,count));
841 : }
842 6 : PetscCall(MatDenseRestoreArray(MS,&S));
843 6 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
844 6 : PetscCall(BVSetActiveColumns(pep->V,0,count));
845 6 : PetscCall(BVDuplicateResize(pep->V,count,&tV));
846 6 : PetscCall(BVCopy(pep->V,tV));
847 : }
848 21 : if (sr->sPres->nconv[0]) {
849 11 : if (divide) {
850 6 : PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[0]));
851 6 : PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[0]));
852 : }
853 94 : for (i=0;i<sr->ndef0;i++) aux[i] = sr->idxDef0[i];
854 61 : for (i=sr->ndef0;i<sr->sPres->nconv[0];i++) aux[i] = sr->indexEig+i-sr->ndef0;
855 11 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
856 11 : PetscCall(MatDenseGetArray(MS,&S));
857 144 : for (i=0;i<sr->sPres->nconv[0];i++) {
858 133 : sr->eigr[aux[i]] = eigr[i];
859 133 : sr->errest[aux[i]] = errest[i];
860 133 : PetscCall(BVGetColumn(pep->V,i,&w));
861 133 : PetscCall(BVInsertVec(sr->V,aux[i],w));
862 133 : PetscCall(BVRestoreColumn(pep->V,i,&w));
863 133 : idx = sr->ld*d*aux[i];
864 133 : PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
865 133 : PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[0]));
866 133 : PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[0]));
867 133 : PetscCall(PetscFree(sr->qinfo[aux[i]].q));
868 133 : PetscCall(PetscMalloc1(sr->sPres->nconv[0],&sr->qinfo[aux[i]].q));
869 133 : PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[0]));
870 133 : sr->qinfo[aux[i]].nq = sr->sPres->nconv[0];
871 : }
872 11 : PetscCall(MatDenseRestoreArray(MS,&S));
873 11 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
874 : }
875 :
876 21 : if (sr->sPres->nconv[1]) {
877 13 : if (divide) {
878 6 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
879 6 : PetscCall(MatDenseGetArray(MS,&S));
880 73 : for (i=0;i<sr->sPres->nconv[1];i++) {
881 67 : PetscCall(PetscArraycpy(S+i*(d*ld),tS+(sr->sPres->nconv[0]+i)*nconv*d,count));
882 67 : PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+(sr->sPres->nconv[0]+i)*nconv*d+nconv,count));
883 : }
884 6 : PetscCall(MatDenseRestoreArray(MS,&S));
885 6 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
886 6 : PetscCall(BVSetActiveColumns(pep->V,0,count));
887 6 : PetscCall(BVCopy(tV,pep->V));
888 6 : PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[1]));
889 6 : PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[1]));
890 : }
891 36 : for (i=0;i<sr->ndef1;i++) aux[i] = sr->idxDef1[i];
892 132 : for (i=sr->ndef1;i<sr->sPres->nconv[1];i++) aux[i] = sr->indexEig+sr->sPres->nconv[0]-sr->ndef0+i-sr->ndef1;
893 13 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
894 13 : PetscCall(MatDenseGetArray(MS,&S));
895 155 : for (i=0;i<sr->sPres->nconv[1];i++) {
896 142 : sr->eigr[aux[i]] = eigr[sr->sPres->nconv[0]+i];
897 142 : sr->errest[aux[i]] = errest[sr->sPres->nconv[0]+i];
898 142 : PetscCall(BVGetColumn(pep->V,i,&w));
899 142 : PetscCall(BVInsertVec(sr->V,aux[i],w));
900 142 : PetscCall(BVRestoreColumn(pep->V,i,&w));
901 142 : idx = sr->ld*d*aux[i];
902 142 : PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
903 142 : PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[1]));
904 142 : PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[1]));
905 142 : PetscCall(PetscFree(sr->qinfo[aux[i]].q));
906 142 : PetscCall(PetscMalloc1(sr->sPres->nconv[1],&sr->qinfo[aux[i]].q));
907 142 : PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[1]));
908 142 : sr->qinfo[aux[i]].nq = sr->sPres->nconv[1];
909 : }
910 13 : PetscCall(MatDenseRestoreArray(MS,&S));
911 13 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
912 : }
913 21 : sPres->neigs = count-sr->ndef0-sr->ndef1;
914 21 : sr->indexEig += sPres->neigs;
915 21 : sPres->nconv[0]-= sr->ndef0;
916 21 : sPres->nconv[1]-= sr->ndef1;
917 21 : PetscCall(PetscFree4(eigr,errest,tS,aux));
918 : } else {
919 1 : sPres->neigs = 0;
920 1 : sPres->nconv[0]= 0;
921 1 : sPres->nconv[1]= 0;
922 : }
923 : /* Global ordering array updating */
924 22 : PetscCall(sortRealEigenvalues(sr->eigr,sr->perm,sr->indexEig,PETSC_FALSE,sr->dir));
925 : /* Check for completion */
926 22 : sPres->comp[0] = PetscNot(sPres->nconv[0] < sPres->nsch[0]);
927 22 : sPres->comp[1] = PetscNot(sPres->nconv[1] < sPres->nsch[1]);
928 22 : PetscCheck(sPres->nconv[0]<=sPres->nsch[0] && sPres->nconv[1]<=sPres->nsch[1],PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
929 22 : if (divide) PetscCall(BVDestroy(&tV));
930 22 : PetscFunctionReturn(PETSC_SUCCESS);
931 : }
932 :
933 22 : static PetscErrorCode PEPLookForDeflation(PEP pep)
934 : {
935 22 : PetscReal val;
936 22 : PetscInt i,count0=0,count1=0;
937 22 : PEP_shift sPres;
938 22 : PetscInt ini,fin;
939 22 : PEP_SR sr;
940 22 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
941 :
942 22 : PetscFunctionBegin;
943 22 : sr = ctx->sr;
944 22 : sPres = sr->sPres;
945 :
946 22 : if (sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
947 : else ini = 0;
948 22 : fin = sr->indexEig;
949 : /* Selection of ends for searching new values */
950 22 : if (!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
951 14 : else sPres->ext[0] = sPres->neighb[0]->value;
952 22 : if (!sPres->neighb[1]) {
953 19 : if (sr->hasEnd) sPres->ext[1] = sr->int1;
954 0 : else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
955 3 : } else sPres->ext[1] = sPres->neighb[1]->value;
956 : /* Selection of values between right and left ends */
957 128 : for (i=ini;i<fin;i++) {
958 108 : val=PetscRealPart(sr->eigr[sr->perm[i]]);
959 : /* Values to the right of left shift */
960 108 : if (sr->dir*(val - sPres->ext[1]) < 0) {
961 106 : if (sr->dir*(val - sPres->value) < 0) count0++;
962 23 : else count1++;
963 : } else break;
964 : }
965 : /* The number of values on each side are found */
966 22 : if (sPres->neighb[0]) {
967 14 : sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
968 14 : PetscCheck(sPres->nsch[0]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
969 8 : } else sPres->nsch[0] = 0;
970 :
971 22 : if (sPres->neighb[1]) {
972 3 : sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
973 3 : PetscCheck(sPres->nsch[1]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
974 19 : } else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);
975 :
976 : /* Completing vector of indexes for deflation */
977 105 : for (i=0;i<count0;i++) sr->idxDef0[i] = sr->perm[ini+i];
978 22 : sr->ndef0 = count0;
979 45 : for (i=0;i<count1;i++) sr->idxDef1[i] = sr->perm[ini+count0+i];
980 22 : sr->ndef1 = count1;
981 22 : PetscFunctionReturn(PETSC_SUCCESS);
982 : }
983 :
984 : /*
985 : Compute a run of Lanczos iterations
986 : */
987 193 : static PetscErrorCode PEPSTOARrun_QSlice(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
988 : {
989 193 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
990 193 : PetscInt i,j,m=*M,l,lock;
991 193 : PetscInt lds,d,ld,offq,nqt,ldds;
992 193 : Vec v=t_[0],t=t_[1],q=t_[2];
993 193 : PetscReal norm,sym=0.0,fro=0.0,*f;
994 193 : PetscScalar *y,*S,sigma;
995 193 : PetscBLASInt j_,one=1;
996 193 : PetscBool lindep;
997 193 : Mat MS;
998 :
999 193 : PetscFunctionBegin;
1000 193 : PetscCall(PetscMalloc1(*M,&y));
1001 193 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
1002 193 : PetscCall(BVTensorGetDegree(ctx->V,&d));
1003 193 : PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
1004 193 : lds = d*ld;
1005 193 : offq = ld;
1006 193 : PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
1007 :
1008 193 : *breakdown = PETSC_FALSE; /* ----- */
1009 193 : PetscCall(STGetShift(pep->st,&sigma));
1010 193 : PetscCall(DSGetDimensions(pep->ds,NULL,&l,NULL,NULL));
1011 193 : PetscCall(BVSetActiveColumns(ctx->V,0,m));
1012 193 : PetscCall(BVSetActiveColumns(pep->V,0,nqt));
1013 3368 : for (j=k;j<m;j++) {
1014 : /* apply operator */
1015 3176 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1016 3176 : PetscCall(MatDenseGetArray(MS,&S));
1017 3176 : PetscCall(BVGetColumn(pep->V,nqt,&t));
1018 3176 : PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+j*lds));
1019 3176 : PetscCall(MatMult(pep->A[1],v,q));
1020 3176 : PetscCall(MatMult(pep->A[2],v,t));
1021 3176 : PetscCall(VecAXPY(q,sigma*pep->sfactor,t));
1022 3176 : PetscCall(VecScale(q,pep->sfactor));
1023 3176 : PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
1024 3176 : PetscCall(MatMult(pep->A[2],v,t));
1025 3176 : PetscCall(VecAXPY(q,pep->sfactor*pep->sfactor,t));
1026 3176 : PetscCall(STMatSolve(pep->st,q,t));
1027 3176 : PetscCall(VecScale(t,-1.0));
1028 3176 : PetscCall(BVRestoreColumn(pep->V,nqt,&t));
1029 :
1030 : /* orthogonalize */
1031 3176 : PetscCall(BVOrthogonalizeColumn(pep->V,nqt,S+(j+1)*lds,&norm,&lindep));
1032 3176 : if (!lindep) {
1033 3176 : *(S+(j+1)*lds+nqt) = norm;
1034 3176 : PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
1035 3176 : nqt++;
1036 : }
1037 88947 : for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
1038 3176 : PetscCall(BVSetActiveColumns(pep->V,0,nqt));
1039 3176 : PetscCall(MatDenseRestoreArray(MS,&S));
1040 3176 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1041 :
1042 : /* level-2 orthogonalization */
1043 3176 : PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
1044 3176 : a[j] = PetscRealPart(y[j]);
1045 3176 : omega[j+1] = (norm > 0)?1.0:-1.0;
1046 3176 : PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
1047 3176 : b[j] = PetscAbsReal(norm);
1048 :
1049 : /* check symmetry */
1050 3176 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&f));
1051 3176 : if (j==k) {
1052 3125 : for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
1053 193 : for (i=0;i<l;i++) y[i] = 0.0;
1054 : }
1055 3176 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&f));
1056 3176 : if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
1057 3176 : PetscCall(PetscBLASIntCast(j,&j_));
1058 3176 : sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
1059 3176 : fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
1060 3176 : if (j>0) fro = SlepcAbs(fro,b[j-1]);
1061 3176 : if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
1062 1 : *symmlost = PETSC_TRUE;
1063 1 : *M=j;
1064 1 : break;
1065 : }
1066 : }
1067 193 : PetscCall(BVSetActiveColumns(pep->V,lock,nqt));
1068 193 : PetscCall(BVSetActiveColumns(ctx->V,0,*M));
1069 193 : PetscCall(PetscFree(y));
1070 193 : PetscFunctionReturn(PETSC_SUCCESS);
1071 : }
1072 :
1073 22 : static PetscErrorCode PEPSTOAR_QSlice(PEP pep,Mat B)
1074 : {
1075 22 : PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
1076 22 : PetscInt j,k,l,nv=0,ld,ldds,t,nq=0,idx;
1077 22 : PetscInt nconv=0,deg=pep->nmat-1,count0=0,count1=0;
1078 22 : PetscScalar *om,sigma,*back,*S,*pQ;
1079 22 : PetscReal beta,norm=1.0,*omega,*a,*b,eta,lambda;
1080 22 : PetscBool breakdown,symmlost=PETSC_FALSE,sinv,falselock=PETSC_TRUE;
1081 22 : Mat MS,MQ,D;
1082 22 : Vec v,vomega;
1083 22 : PEP_SR sr;
1084 22 : BVOrthogType otype;
1085 22 : BVOrthogBlockType obtype;
1086 :
1087 22 : PetscFunctionBegin;
1088 : /* Resize if needed for deflating vectors */
1089 22 : sr = ctx->sr;
1090 22 : sigma = sr->sPres->value;
1091 22 : k = sr->ndef0+sr->ndef1;
1092 22 : pep->ncv = ctx->ncv+k;
1093 22 : pep->nev = ctx->nev+k;
1094 22 : PetscCall(PEPAllocateSolution(pep,3));
1095 22 : PetscCall(BVDestroy(&ctx->V));
1096 22 : PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
1097 22 : PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
1098 22 : PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
1099 22 : PetscCall(DSAllocate(pep->ds,pep->ncv+2));
1100 22 : PetscCall(PetscMalloc1(pep->ncv,&back));
1101 22 : PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
1102 22 : PetscCall(BVSetMatrix(ctx->V,B,PETSC_TRUE));
1103 22 : PetscCheck(ctx->lock,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"A locking variant is needed for spectrum slicing");
1104 : /* undocumented option to use a cheaper locking instead of the true locking */
1105 22 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL));
1106 22 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
1107 22 : PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
1108 22 : PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
1109 :
1110 : /* Get the starting Arnoldi vector */
1111 22 : PetscCall(BVSetActiveColumns(pep->V,0,1));
1112 22 : PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
1113 22 : PetscCall(BVSetActiveColumns(ctx->V,0,1));
1114 22 : if (k) {
1115 : /* Insert deflated vectors */
1116 10 : PetscCall(BVSetActiveColumns(pep->V,0,0));
1117 10 : idx = sr->ndef0?sr->idxDef0[0]:sr->idxDef1[0];
1118 116 : for (j=0;j<k;j++) {
1119 106 : PetscCall(BVGetColumn(pep->V,j,&v));
1120 106 : PetscCall(BVCopyVec(sr->V,sr->qinfo[idx].q[j],v));
1121 106 : PetscCall(BVRestoreColumn(pep->V,j,&v));
1122 : }
1123 : /* Update innerproduct matrix */
1124 10 : PetscCall(BVSetActiveColumns(ctx->V,0,0));
1125 10 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1126 10 : PetscCall(BVSetActiveColumns(pep->V,0,k));
1127 10 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1128 :
1129 10 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
1130 10 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1131 10 : PetscCall(MatDenseGetArray(MS,&S));
1132 93 : for (j=0;j<sr->ndef0;j++) {
1133 83 : PetscCall(PetscArrayzero(S+j*ld*deg,ld*deg));
1134 83 : PetscCall(PetscArraycpy(S+j*ld*deg,sr->S+sr->idxDef0[j]*sr->ld*deg,k));
1135 83 : PetscCall(PetscArraycpy(S+j*ld*deg+ld,sr->S+sr->idxDef0[j]*sr->ld*deg+sr->ld,k));
1136 83 : pep->eigr[j] = sr->eigr[sr->idxDef0[j]];
1137 83 : pep->errest[j] = sr->errest[sr->idxDef0[j]];
1138 : }
1139 33 : for (j=0;j<sr->ndef1;j++) {
1140 23 : PetscCall(PetscArrayzero(S+(j+sr->ndef0)*ld*deg,ld*deg));
1141 23 : PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg,sr->S+sr->idxDef1[j]*sr->ld*deg,k));
1142 23 : PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg+ld,sr->S+sr->idxDef1[j]*sr->ld*deg+sr->ld,k));
1143 23 : pep->eigr[j+sr->ndef0] = sr->eigr[sr->idxDef1[j]];
1144 23 : pep->errest[j+sr->ndef0] = sr->errest[sr->idxDef1[j]];
1145 : }
1146 10 : PetscCall(MatDenseRestoreArray(MS,&S));
1147 10 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1148 10 : PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
1149 10 : PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+1,&vomega));
1150 10 : PetscCall(VecGetArray(vomega,&om));
1151 116 : for (j=0;j<k;j++) {
1152 106 : PetscCall(BVOrthogonalizeColumn(ctx->V,j,NULL,&norm,NULL));
1153 106 : PetscCall(BVScaleColumn(ctx->V,j,1/norm));
1154 125 : om[j] = (norm>=0.0)?1.0:-1.0;
1155 : }
1156 10 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1157 10 : PetscCall(MatDenseGetArray(MS,&S));
1158 30 : for (j=0;j<deg;j++) {
1159 20 : PetscCall(BVSetRandomColumn(pep->V,k+j));
1160 20 : PetscCall(BVOrthogonalizeColumn(pep->V,k+j,S+k*ld*deg+j*ld,&norm,NULL));
1161 20 : PetscCall(BVScaleColumn(pep->V,k+j,1.0/norm));
1162 20 : S[k*ld*deg+j*ld+k+j] = norm;
1163 : }
1164 10 : PetscCall(MatDenseRestoreArray(MS,&S));
1165 10 : PetscCall(BVSetActiveColumns(pep->V,0,k+deg));
1166 10 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1167 10 : PetscCall(BVOrthogonalizeColumn(ctx->V,k,NULL,&norm,NULL));
1168 10 : PetscCall(BVScaleColumn(ctx->V,k,1.0/norm));
1169 10 : om[k] = (norm>=0.0)?1.0:-1.0;
1170 10 : PetscCall(VecRestoreArray(vomega,&om));
1171 10 : PetscCall(BVSetSignature(ctx->V,vomega));
1172 10 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
1173 10 : PetscCall(VecGetArray(vomega,&om));
1174 116 : for (j=0;j<k;j++) a[j] = PetscRealPart(om[j]/(pep->eigr[j]-sigma));
1175 10 : PetscCall(VecRestoreArray(vomega,&om));
1176 10 : PetscCall(VecDestroy(&vomega));
1177 10 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
1178 10 : PetscCall(DSGetArray(pep->ds,DS_MAT_Q,&pQ));
1179 10 : PetscCall(PetscArrayzero(pQ,ldds*k));
1180 116 : for (j=0;j<k;j++) pQ[j+j*ldds] = 1.0;
1181 10 : PetscCall(DSRestoreArray(pep->ds,DS_MAT_Q,&pQ));
1182 : }
1183 22 : PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
1184 22 : PetscCall(DSSetDimensions(pep->ds,k+1,PETSC_DETERMINE,PETSC_DETERMINE));
1185 22 : PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1186 22 : PetscCall(BVGetSignature(ctx->V,vomega));
1187 22 : PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1188 :
1189 22 : PetscCall(PetscInfo(pep,"Start STOAR: sigma=%g in [%g,%g], for deflation: left=%" PetscInt_FMT " right=%" PetscInt_FMT ", searching: left=%" PetscInt_FMT " right=%" PetscInt_FMT "\n",(double)sr->sPres->value,(double)(sr->sPres->neighb[0]?sr->sPres->neighb[0]->value:sr->int0),(double)(sr->sPres->neighb[1]?sr->sPres->neighb[1]->value:sr->int1),sr->ndef0,sr->ndef1,sr->sPres->nsch[0],sr->sPres->nsch[1]));
1190 :
1191 : /* Restart loop */
1192 22 : l = 0;
1193 22 : pep->nconv = k;
1194 22 : while (pep->reason == PEP_CONVERGED_ITERATING) {
1195 193 : pep->its++;
1196 193 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
1197 193 : b = a+ldds;
1198 193 : PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));
1199 :
1200 : /* Compute an nv-step Lanczos factorization */
1201 193 : nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
1202 193 : PetscCall(PEPSTOARrun_QSlice(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work));
1203 193 : beta = b[nv-1];
1204 193 : if (symmlost && nv==pep->nconv+l) {
1205 1 : pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
1206 1 : pep->nconv = nconv;
1207 1 : PetscCall(PetscInfo(pep,"Symmetry lost in STOAR sigma=%g nconv=%" PetscInt_FMT "\n",(double)sr->sPres->value,nconv));
1208 1 : if (falselock || !ctx->lock) {
1209 1 : PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
1210 1 : PetscCall(BVTensorCompress(ctx->V,0));
1211 : }
1212 : break;
1213 : }
1214 192 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
1215 192 : PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega));
1216 192 : PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
1217 192 : if (l==0) PetscCall(DSSetState(pep->ds,DS_STATE_INTERMEDIATE));
1218 170 : else PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
1219 :
1220 : /* Solve projected problem */
1221 192 : PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
1222 192 : PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
1223 192 : PetscCall(DSUpdateExtraRow(pep->ds));
1224 192 : PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
1225 :
1226 : /* Check convergence */
1227 : /* PetscCall(PEPSTOARpreKConvergence(pep,nv,&norm,pep->work));*/
1228 192 : norm = 1.0;
1229 192 : PetscCall(DSGetDimensions(pep->ds,NULL,NULL,NULL,&t));
1230 192 : PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k));
1231 192 : PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
1232 1231 : for (j=0;j<k;j++) back[j] = pep->eigr[j];
1233 192 : PetscCall(STBackTransform(pep->st,k,back,pep->eigi));
1234 : count0=count1=0;
1235 1231 : for (j=0;j<k;j++) {
1236 1039 : lambda = PetscRealPart(back[j]);
1237 1039 : if ((sr->dir*(sr->sPres->value - lambda) > 0) && (sr->dir*(lambda - sr->sPres->ext[0]) > 0)) count0++;
1238 1039 : if ((sr->dir*(lambda - sr->sPres->value) > 0) && (sr->dir*(sr->sPres->ext[1] - lambda) > 0)) count1++;
1239 : }
1240 192 : if ((count0-sr->ndef0 >= sr->sPres->nsch[0]) && (count1-sr->ndef1 >= sr->sPres->nsch[1])) pep->reason = PEP_CONVERGED_TOL;
1241 : /* Update l */
1242 192 : if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
1243 : else {
1244 171 : l = PetscMax(1,(PetscInt)((nv-k)/2));
1245 171 : l = PetscMin(l,t);
1246 171 : PetscCall(DSGetTruncateSize(pep->ds,k,t,&l));
1247 171 : if (!breakdown) {
1248 : /* Prepare the Rayleigh quotient for restart */
1249 171 : PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
1250 : }
1251 : }
1252 192 : nconv = k;
1253 192 : if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
1254 192 : if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1255 :
1256 : /* Update S */
1257 192 : PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
1258 192 : PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
1259 192 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
1260 :
1261 : /* Copy last column of S */
1262 192 : PetscCall(BVCopyColumn(ctx->V,nv,k+l));
1263 192 : PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
1264 192 : if (k+l) {
1265 191 : PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DETERMINE,PETSC_DETERMINE));
1266 191 : PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1267 191 : PetscCall(BVSetSignature(ctx->V,vomega));
1268 191 : PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1269 : }
1270 :
1271 192 : if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
1272 : /* stop if breakdown */
1273 0 : PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
1274 0 : pep->reason = PEP_DIVERGED_BREAKDOWN;
1275 : }
1276 192 : if (pep->reason != PEP_CONVERGED_ITERATING) l--;
1277 192 : PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
1278 192 : if (k+l+deg<=nq) {
1279 192 : PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
1280 192 : if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
1281 192 : else PetscCall(BVTensorCompress(ctx->V,0));
1282 : }
1283 192 : pep->nconv = k;
1284 214 : PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
1285 : }
1286 22 : sr->itsKs += pep->its;
1287 22 : if (pep->nconv>0) {
1288 21 : PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
1289 21 : PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
1290 21 : PetscCall(BVSetActiveColumns(pep->V,0,nq));
1291 21 : if (nq>pep->nconv) {
1292 21 : PetscCall(BVTensorCompress(ctx->V,pep->nconv));
1293 21 : PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
1294 : }
1295 450 : for (j=0;j<pep->nconv;j++) {
1296 429 : pep->eigr[j] *= pep->sfactor;
1297 429 : pep->eigi[j] *= pep->sfactor;
1298 : }
1299 : }
1300 22 : PetscCall(PetscInfo(pep,"Finished STOAR: nconv=%" PetscInt_FMT " (deflated=%" PetscInt_FMT ", left=%" PetscInt_FMT ", right=%" PetscInt_FMT ")\n",pep->nconv,sr->ndef0+sr->ndef1,count0-sr->ndef0,count1-sr->ndef1));
1301 22 : PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
1302 22 : PetscCall(RGPopScale(pep->rg));
1303 :
1304 22 : PetscCheck(pep->reason!=PEP_DIVERGED_SYMMETRY_LOST || nconv>=sr->ndef0+sr->ndef1,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
1305 22 : if (pep->reason == PEP_DIVERGED_SYMMETRY_LOST && nconv==sr->ndef0+sr->ndef1) {
1306 0 : PetscCheck(++sr->symmlost<=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
1307 22 : } else sr->symmlost = 0;
1308 :
1309 22 : PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
1310 22 : PetscCall(PetscFree(back));
1311 22 : PetscFunctionReturn(PETSC_SUCCESS);
1312 : }
1313 :
1314 7 : static PetscErrorCode PEPQSliceGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
1315 : {
1316 7 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
1317 7 : PEP_SR sr=ctx->sr;
1318 7 : PetscInt i=0,j,tmpi;
1319 7 : PetscReal v,tmpr;
1320 7 : PEP_shift s;
1321 :
1322 7 : PetscFunctionBegin;
1323 7 : PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
1324 7 : PetscCheck(sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
1325 7 : if (!sr->s0) { /* PEPSolve not called yet */
1326 0 : *n = 2;
1327 : } else {
1328 7 : *n = 1;
1329 7 : s = sr->s0;
1330 28 : while (s) {
1331 21 : (*n)++;
1332 21 : s = s->neighb[1];
1333 : }
1334 : }
1335 7 : PetscCall(PetscMalloc1(*n,shifts));
1336 7 : PetscCall(PetscMalloc1(*n,inertias));
1337 7 : if (!sr->s0) { /* PEPSolve not called yet */
1338 0 : (*shifts)[0] = sr->int0;
1339 0 : (*shifts)[1] = sr->int1;
1340 0 : (*inertias)[0] = sr->inertia0;
1341 0 : (*inertias)[1] = sr->inertia1;
1342 : } else {
1343 : s = sr->s0;
1344 28 : while (s) {
1345 21 : (*shifts)[i] = s->value;
1346 21 : (*inertias)[i++] = s->inertia;
1347 21 : s = s->neighb[1];
1348 : }
1349 7 : (*shifts)[i] = sr->int1;
1350 7 : (*inertias)[i] = sr->inertia1;
1351 : }
1352 : /* remove possible duplicate in last position */
1353 7 : if ((*shifts)[(*n)-1]==(*shifts)[(*n)-2]) (*n)--;
1354 : /* sort result */
1355 32 : for (i=0;i<*n;i++) {
1356 25 : v = (*shifts)[i];
1357 81 : for (j=i+1;j<*n;j++) {
1358 56 : if (v > (*shifts)[j]) {
1359 20 : SlepcSwap((*shifts)[i],(*shifts)[j],tmpr);
1360 20 : SlepcSwap((*inertias)[i],(*inertias)[j],tmpi);
1361 20 : v = (*shifts)[i];
1362 : }
1363 : }
1364 : }
1365 7 : PetscFunctionReturn(PETSC_SUCCESS);
1366 : }
1367 :
1368 8 : PetscErrorCode PEPSolve_STOAR_QSlice(PEP pep)
1369 : {
1370 8 : PetscInt i,j,ti,deg=pep->nmat-1;
1371 8 : PetscReal newS;
1372 8 : PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
1373 8 : PEP_SR sr=ctx->sr;
1374 8 : Mat S,B;
1375 8 : PetscScalar *pS;
1376 :
1377 8 : PetscFunctionBegin;
1378 8 : PetscCall(PetscCitationsRegister(citation,&cited));
1379 :
1380 : /* Only with eigenvalues present in the interval ...*/
1381 8 : if (sr->numEigs==0) {
1382 1 : pep->reason = PEP_CONVERGED_TOL;
1383 1 : PetscFunctionReturn(PETSC_SUCCESS);
1384 : }
1385 :
1386 : /* Inner product matrix */
1387 7 : PetscCall(PEPSTOARSetUpInnerMatrix(pep,&B));
1388 :
1389 : /* Array of pending shifts */
1390 7 : sr->maxPend = 100; /* Initial size */
1391 7 : sr->nPend = 0;
1392 7 : PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
1393 7 : PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
1394 : /* extract first shift */
1395 7 : sr->sPrev = NULL;
1396 7 : sr->sPres = sr->pending[--sr->nPend];
1397 7 : sr->sPres->inertia = sr->inertia0;
1398 7 : pep->target = sr->sPres->value;
1399 7 : sr->s0 = sr->sPres;
1400 7 : sr->indexEig = 0;
1401 :
1402 176 : for (i=0;i<sr->numEigs;i++) {
1403 169 : sr->eigr[i] = 0.0;
1404 169 : sr->eigi[i] = 0.0;
1405 169 : sr->errest[i] = 0.0;
1406 169 : sr->perm[i] = i;
1407 : }
1408 : /* Vectors for deflation */
1409 7 : PetscCall(PetscMalloc2(sr->numEigs,&sr->idxDef0,sr->numEigs,&sr->idxDef1));
1410 7 : sr->indexEig = 0;
1411 7 : while (sr->sPres) {
1412 : /* Search for deflation */
1413 22 : PetscCall(PEPLookForDeflation(pep));
1414 : /* KrylovSchur */
1415 22 : PetscCall(PEPSTOAR_QSlice(pep,B));
1416 :
1417 22 : PetscCall(PEPStoreEigenpairs(pep));
1418 : /* Select new shift */
1419 22 : if (!sr->sPres->comp[1]) {
1420 12 : PetscCall(PEPGetNewShiftValue(pep,1,&newS));
1421 12 : PetscCall(PEPCreateShift(pep,newS,sr->sPres,sr->sPres->neighb[1]));
1422 : }
1423 22 : if (!sr->sPres->comp[0]) {
1424 : /* Completing earlier interval */
1425 3 : PetscCall(PEPGetNewShiftValue(pep,0,&newS));
1426 3 : PetscCall(PEPCreateShift(pep,newS,sr->sPres->neighb[0],sr->sPres));
1427 : }
1428 : /* Preparing for a new search of values */
1429 29 : PetscCall(PEPExtractShift(pep));
1430 : }
1431 :
1432 : /* Updating pep values prior to exit */
1433 7 : PetscCall(PetscFree2(sr->idxDef0,sr->idxDef1));
1434 7 : PetscCall(PetscFree(sr->pending));
1435 7 : pep->nconv = sr->indexEig;
1436 7 : pep->reason = PEP_CONVERGED_TOL;
1437 7 : pep->its = sr->itsKs;
1438 7 : pep->nev = sr->indexEig;
1439 7 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,pep->nconv,pep->nconv,NULL,&S));
1440 7 : PetscCall(MatDenseGetArray(S,&pS));
1441 176 : for (i=0;i<pep->nconv;i++) {
1442 2466 : for (j=0;j<sr->qinfo[i].nq;j++) pS[i*pep->nconv+sr->qinfo[i].q[j]] = *(sr->S+i*sr->ld*deg+j);
1443 : }
1444 7 : PetscCall(MatDenseRestoreArray(S,&pS));
1445 7 : PetscCall(BVSetActiveColumns(sr->V,0,pep->nconv));
1446 7 : PetscCall(BVMultInPlace(sr->V,S,0,pep->nconv));
1447 7 : PetscCall(MatDestroy(&S));
1448 7 : PetscCall(BVDestroy(&pep->V));
1449 7 : pep->V = sr->V;
1450 7 : PetscCall(PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm));
1451 7 : pep->eigr = sr->eigr;
1452 7 : pep->eigi = sr->eigi;
1453 7 : pep->perm = sr->perm;
1454 7 : pep->errest = sr->errest;
1455 7 : if (sr->dir<0) {
1456 48 : for (i=0;i<pep->nconv/2;i++) {
1457 42 : ti = sr->perm[i]; sr->perm[i] = sr->perm[pep->nconv-1-i]; sr->perm[pep->nconv-1-i] = ti;
1458 : }
1459 : }
1460 7 : PetscCall(PetscFree(ctx->inertias));
1461 7 : PetscCall(PetscFree(ctx->shifts));
1462 7 : PetscCall(MatDestroy(&B));
1463 7 : PetscCall(PEPQSliceGetInertias(pep,&ctx->nshifts,&ctx->shifts,&ctx->inertias));
1464 7 : PetscFunctionReturn(PETSC_SUCCESS);
1465 : }
|
