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 quadratic eigensolver: "qarnoldi"
12 :
13 : Method: Q-Arnoldi
14 :
15 : Algorithm:
16 :
17 : Quadratic Arnoldi with Krylov-Schur type restart.
18 :
19 : References:
20 :
21 : [1] K. Meerbergen, "The Quadratic Arnoldi method for the solution
22 : of the quadratic eigenvalue problem", SIAM J. Matrix Anal.
23 : Appl. 30(4):1462-1482, 2008.
24 : */
25 :
26 : #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
27 : #include <petscblaslapack.h>
28 :
29 : typedef struct {
30 : PetscReal keep; /* restart parameter */
31 : PetscBool lock; /* locking/non-locking variant */
32 : } PEP_QARNOLDI;
33 :
34 17 : static PetscErrorCode PEPSetUp_QArnoldi(PEP pep)
35 : {
36 17 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
37 17 : PetscBool flg;
38 :
39 17 : PetscFunctionBegin;
40 17 : PEPCheckQuadratic(pep);
41 17 : PEPCheckShiftSinvert(pep);
42 17 : PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
43 17 : PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
44 17 : if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);
45 17 : if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
46 17 : PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
47 :
48 17 : PetscCall(STGetTransform(pep->st,&flg));
49 17 : PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");
50 :
51 : /* set default extraction */
52 17 : if (!pep->extract) pep->extract = PEP_EXTRACT_NONE;
53 17 : PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_EXTRACT);
54 :
55 17 : if (!ctx->keep) ctx->keep = 0.5;
56 :
57 17 : PetscCall(PEPAllocateSolution(pep,0));
58 17 : PetscCall(PEPSetWorkVecs(pep,4));
59 :
60 17 : PetscCall(DSSetType(pep->ds,DSNHEP));
61 17 : PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
62 17 : PetscCall(DSAllocate(pep->ds,pep->ncv+1));
63 17 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 17 : static PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
67 : {
68 17 : PetscInt k=pep->nconv;
69 17 : Mat X,X0;
70 :
71 17 : PetscFunctionBegin;
72 17 : if (pep->nconv==0) PetscFunctionReturn(PETSC_SUCCESS);
73 17 : PetscCall(DSVectors(pep->ds,DS_MAT_X,NULL,NULL));
74 :
75 : /* update vectors V = V*X */
76 17 : PetscCall(DSGetMat(pep->ds,DS_MAT_X,&X));
77 17 : PetscCall(MatDenseGetSubMatrix(X,0,k,0,k,&X0));
78 17 : PetscCall(BVMultInPlace(pep->V,X0,0,k));
79 17 : PetscCall(MatDenseRestoreSubMatrix(X,&X0));
80 17 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_X,&X));
81 17 : PetscFunctionReturn(PETSC_SUCCESS);
82 : }
83 :
84 : /*
85 : Compute a step of Classical Gram-Schmidt orthogonalization
86 : */
87 3109 : static PetscErrorCode PEPQArnoldiCGS(PEP pep,PetscScalar *H,PetscBLASInt ldh,PetscScalar *h,PetscBLASInt j,BV V,Vec t,Vec v,Vec w,PetscReal *onorm,PetscReal *norm,PetscScalar *work)
88 : {
89 3109 : PetscBLASInt ione = 1,j_1 = j+1;
90 3109 : PetscReal x,y;
91 3109 : PetscScalar dot,one = 1.0,zero = 0.0;
92 :
93 3109 : PetscFunctionBegin;
94 : /* compute norm of v and w */
95 3109 : if (onorm) {
96 1560 : PetscCall(VecNorm(v,NORM_2,&x));
97 1560 : PetscCall(VecNorm(w,NORM_2,&y));
98 1560 : *onorm = SlepcAbs(x,y);
99 : }
100 :
101 : /* orthogonalize: compute h */
102 3109 : PetscCall(BVDotVec(V,v,h));
103 3109 : PetscCall(BVDotVec(V,w,work));
104 3109 : if (j>0) PetscCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
105 3109 : PetscCall(VecDot(w,t,&dot));
106 3109 : h[j] += dot;
107 :
108 : /* orthogonalize: update v and w */
109 3109 : PetscCall(BVMultVec(V,-1.0,1.0,v,h));
110 3109 : if (j>0) {
111 3090 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
112 3090 : PetscCall(BVMultVec(V,-1.0,1.0,w,work));
113 : }
114 3109 : PetscCall(VecAXPY(w,-h[j],t));
115 :
116 : /* compute norm of v and w */
117 3109 : if (norm) {
118 2979 : PetscCall(VecNorm(v,NORM_2,&x));
119 2979 : PetscCall(VecNorm(w,NORM_2,&y));
120 2979 : *norm = SlepcAbs(x,y);
121 : }
122 3109 : PetscFunctionReturn(PETSC_SUCCESS);
123 : }
124 :
125 : /*
126 : Compute a run of Q-Arnoldi iterations
127 : */
128 128 : static PetscErrorCode PEPQArnoldi(PEP pep,Mat A,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
129 : {
130 128 : PetscInt i,j,l,m = *M,ldh;
131 128 : PetscBLASInt jj,ldhh;
132 128 : Vec t = pep->work[2],u = pep->work[3];
133 128 : BVOrthogRefineType refinement;
134 128 : PetscReal norm=0.0,onorm,eta;
135 128 : PetscScalar *H,*c = work + m;
136 :
137 128 : PetscFunctionBegin;
138 128 : *beta = 0.0;
139 128 : PetscCall(MatDenseGetArray(A,&H));
140 128 : PetscCall(MatDenseGetLDA(A,&ldh));
141 128 : PetscCall(BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL));
142 128 : PetscCall(BVInsertVec(pep->V,k,v));
143 1747 : for (j=k;j<m;j++) {
144 : /* apply operator */
145 1619 : PetscCall(VecCopy(w,t));
146 1619 : if (pep->Dr) PetscCall(VecPointwiseMult(v,v,pep->Dr));
147 1619 : PetscCall(STMatMult(pep->st,0,v,u));
148 1619 : PetscCall(VecCopy(t,v));
149 1619 : if (pep->Dr) PetscCall(VecPointwiseMult(t,t,pep->Dr));
150 1619 : PetscCall(STMatMult(pep->st,1,t,w));
151 1619 : PetscCall(VecAXPY(u,pep->sfactor,w));
152 1619 : PetscCall(STMatSolve(pep->st,u,w));
153 1619 : PetscCall(VecScale(w,-1.0/(pep->sfactor*pep->sfactor)));
154 1619 : if (pep->Dr) PetscCall(VecPointwiseDivide(w,w,pep->Dr));
155 1619 : PetscCall(VecCopy(v,t));
156 1619 : PetscCall(BVSetActiveColumns(pep->V,0,j+1));
157 :
158 : /* orthogonalize */
159 1619 : PetscCall(PetscBLASIntCast(j,&jj));
160 1619 : PetscCall(PetscBLASIntCast(ldh,&ldhh));
161 1619 : switch (refinement) {
162 59 : case BV_ORTHOG_REFINE_NEVER:
163 59 : PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,&norm,work));
164 59 : *breakdown = PETSC_FALSE;
165 59 : break;
166 130 : case BV_ORTHOG_REFINE_ALWAYS:
167 130 : PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,NULL,work));
168 130 : PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,&onorm,&norm,work));
169 2045 : for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
170 130 : if (norm < eta * onorm) *breakdown = PETSC_TRUE;
171 130 : else *breakdown = PETSC_FALSE;
172 : break;
173 1430 : case BV_ORTHOG_REFINE_IFNEEDED:
174 1430 : PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,&onorm,&norm,work));
175 : /* ||q|| < eta ||h|| */
176 : l = 1;
177 2790 : while (l<3 && norm < eta * onorm) {
178 1360 : l++;
179 1360 : onorm = norm;
180 1360 : PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,NULL,&norm,work));
181 24853 : for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
182 : }
183 1430 : if (norm < eta * onorm) *breakdown = PETSC_TRUE;
184 1430 : else *breakdown = PETSC_FALSE;
185 : break;
186 : }
187 1619 : PetscCall(VecScale(v,1.0/norm));
188 1619 : PetscCall(VecScale(w,1.0/norm));
189 :
190 1619 : H[j+1+ldh*j] = norm;
191 1619 : if (j<m-1) PetscCall(BVInsertVec(pep->V,j+1,v));
192 : }
193 128 : *beta = norm;
194 128 : PetscCall(MatDenseRestoreArray(A,&H));
195 128 : PetscFunctionReturn(PETSC_SUCCESS);
196 : }
197 :
198 17 : static PetscErrorCode PEPSolve_QArnoldi(PEP pep)
199 : {
200 17 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
201 17 : PetscInt j,k,l,lwork,nv,nconv;
202 17 : Vec v=pep->work[0],w=pep->work[1];
203 17 : Mat Q,S;
204 17 : PetscScalar *work;
205 17 : PetscReal beta,norm,x,y;
206 17 : PetscBool breakdown=PETSC_FALSE,sinv;
207 :
208 17 : PetscFunctionBegin;
209 17 : lwork = 7*pep->ncv;
210 17 : PetscCall(PetscMalloc1(lwork,&work));
211 17 : PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
212 17 : PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
213 17 : PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
214 :
215 : /* Get the starting Arnoldi vector */
216 51 : for (j=0;j<2;j++) {
217 34 : if (j>=pep->nini) PetscCall(BVSetRandomColumn(pep->V,j));
218 : }
219 17 : PetscCall(BVCopyVec(pep->V,0,v));
220 17 : PetscCall(BVCopyVec(pep->V,1,w));
221 17 : PetscCall(VecNorm(v,NORM_2,&x));
222 17 : PetscCall(VecNorm(w,NORM_2,&y));
223 17 : norm = SlepcAbs(x,y);
224 17 : PetscCall(VecScale(v,1.0/norm));
225 17 : PetscCall(VecScale(w,1.0/norm));
226 :
227 : /* clean projected matrix (including the extra-arrow) */
228 17 : PetscCall(DSSetDimensions(pep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
229 17 : PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
230 17 : PetscCall(MatZeroEntries(S));
231 17 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));
232 :
233 : /* Restart loop */
234 17 : l = 0;
235 17 : while (pep->reason == PEP_CONVERGED_ITERATING) {
236 128 : pep->its++;
237 :
238 : /* Compute an nv-step Arnoldi factorization */
239 128 : nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
240 128 : PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
241 128 : PetscCall(PEPQArnoldi(pep,S,pep->nconv+l,&nv,v,w,&beta,&breakdown,work));
242 128 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));
243 128 : PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
244 128 : PetscCall(DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
245 128 : PetscCall(BVSetActiveColumns(pep->V,pep->nconv,nv));
246 :
247 : /* Solve projected problem */
248 128 : PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
249 128 : PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
250 128 : PetscCall(DSUpdateExtraRow(pep->ds));
251 128 : PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
252 :
253 : /* Check convergence */
254 128 : PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k));
255 128 : PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
256 128 : nconv = k;
257 :
258 : /* Update l */
259 128 : if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
260 : else {
261 111 : l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
262 111 : PetscCall(DSGetTruncateSize(pep->ds,k,nv,&l));
263 : }
264 128 : if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
265 128 : if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
266 :
267 128 : if (pep->reason == PEP_CONVERGED_ITERATING) {
268 111 : if (PetscUnlikely(breakdown)) {
269 : /* Stop if breakdown */
270 0 : PetscCall(PetscInfo(pep,"Breakdown Quadratic Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
271 0 : pep->reason = PEP_DIVERGED_BREAKDOWN;
272 : } else {
273 : /* Prepare the Rayleigh quotient for restart */
274 111 : PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
275 : }
276 : }
277 : /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
278 128 : PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&Q));
279 128 : PetscCall(BVMultInPlace(pep->V,Q,pep->nconv,k+l));
280 128 : PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&Q));
281 :
282 128 : pep->nconv = k;
283 145 : PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
284 : }
285 17 : PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
286 98 : for (j=0;j<pep->nconv;j++) {
287 81 : pep->eigr[j] *= pep->sfactor;
288 81 : pep->eigi[j] *= pep->sfactor;
289 : }
290 :
291 17 : PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
292 17 : PetscCall(RGPopScale(pep->rg));
293 :
294 17 : PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
295 17 : PetscCall(PetscFree(work));
296 17 : PetscFunctionReturn(PETSC_SUCCESS);
297 : }
298 :
299 1 : static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
300 : {
301 1 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
302 :
303 1 : PetscFunctionBegin;
304 1 : if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
305 : else {
306 1 : PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
307 1 : ctx->keep = keep;
308 : }
309 1 : PetscFunctionReturn(PETSC_SUCCESS);
310 : }
311 :
312 : /*@
313 : PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
314 : method, in particular the proportion of basis vectors that must be kept
315 : after restart.
316 :
317 : Logically Collective
318 :
319 : Input Parameters:
320 : + pep - the eigenproblem solver context
321 : - keep - the number of vectors to be kept at restart
322 :
323 : Options Database Key:
324 : . -pep_qarnoldi_restart - Sets the restart parameter
325 :
326 : Notes:
327 : Allowed values are in the range [0.1,0.9]. The default is 0.5.
328 :
329 : Level: advanced
330 :
331 : .seealso: PEPQArnoldiGetRestart()
332 : @*/
333 1 : PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
334 : {
335 1 : PetscFunctionBegin;
336 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
337 3 : PetscValidLogicalCollectiveReal(pep,keep,2);
338 1 : PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
339 1 : PetscFunctionReturn(PETSC_SUCCESS);
340 : }
341 :
342 1 : static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
343 : {
344 1 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
345 :
346 1 : PetscFunctionBegin;
347 1 : *keep = ctx->keep;
348 1 : PetscFunctionReturn(PETSC_SUCCESS);
349 : }
350 :
351 : /*@
352 : PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.
353 :
354 : Not Collective
355 :
356 : Input Parameter:
357 : . pep - the eigenproblem solver context
358 :
359 : Output Parameter:
360 : . keep - the restart parameter
361 :
362 : Level: advanced
363 :
364 : .seealso: PEPQArnoldiSetRestart()
365 : @*/
366 1 : PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
367 : {
368 1 : PetscFunctionBegin;
369 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
370 1 : PetscAssertPointer(keep,2);
371 1 : PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
372 1 : PetscFunctionReturn(PETSC_SUCCESS);
373 : }
374 :
375 2 : static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
376 : {
377 2 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
378 :
379 2 : PetscFunctionBegin;
380 2 : ctx->lock = lock;
381 2 : PetscFunctionReturn(PETSC_SUCCESS);
382 : }
383 :
384 : /*@
385 : PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
386 : the Q-Arnoldi method.
387 :
388 : Logically Collective
389 :
390 : Input Parameters:
391 : + pep - the eigenproblem solver context
392 : - lock - true if the locking variant must be selected
393 :
394 : Options Database Key:
395 : . -pep_qarnoldi_locking - Sets the locking flag
396 :
397 : Notes:
398 : The default is to lock converged eigenpairs when the method restarts.
399 : This behaviour can be changed so that all directions are kept in the
400 : working subspace even if already converged to working accuracy (the
401 : non-locking variant).
402 :
403 : Level: advanced
404 :
405 : .seealso: PEPQArnoldiGetLocking()
406 : @*/
407 2 : PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
408 : {
409 2 : PetscFunctionBegin;
410 2 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
411 6 : PetscValidLogicalCollectiveBool(pep,lock,2);
412 2 : PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
413 2 : PetscFunctionReturn(PETSC_SUCCESS);
414 : }
415 :
416 1 : static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
417 : {
418 1 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
419 :
420 1 : PetscFunctionBegin;
421 1 : *lock = ctx->lock;
422 1 : PetscFunctionReturn(PETSC_SUCCESS);
423 : }
424 :
425 : /*@
426 : PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.
427 :
428 : Not Collective
429 :
430 : Input Parameter:
431 : . pep - the eigenproblem solver context
432 :
433 : Output Parameter:
434 : . lock - the locking flag
435 :
436 : Level: advanced
437 :
438 : .seealso: PEPQArnoldiSetLocking()
439 : @*/
440 1 : PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
441 : {
442 1 : PetscFunctionBegin;
443 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
444 1 : PetscAssertPointer(lock,2);
445 1 : PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
446 1 : PetscFunctionReturn(PETSC_SUCCESS);
447 : }
448 :
449 15 : static PetscErrorCode PEPSetFromOptions_QArnoldi(PEP pep,PetscOptionItems *PetscOptionsObject)
450 : {
451 15 : PetscBool flg,lock;
452 15 : PetscReal keep;
453 :
454 15 : PetscFunctionBegin;
455 15 : PetscOptionsHeadBegin(PetscOptionsObject,"PEP Q-Arnoldi Options");
456 :
457 15 : PetscCall(PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg));
458 15 : if (flg) PetscCall(PEPQArnoldiSetRestart(pep,keep));
459 :
460 15 : PetscCall(PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg));
461 15 : if (flg) PetscCall(PEPQArnoldiSetLocking(pep,lock));
462 :
463 15 : PetscOptionsHeadEnd();
464 15 : PetscFunctionReturn(PETSC_SUCCESS);
465 : }
466 :
467 0 : static PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
468 : {
469 0 : PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
470 0 : PetscBool isascii;
471 :
472 0 : PetscFunctionBegin;
473 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
474 0 : if (isascii) {
475 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
476 0 : PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
477 : }
478 0 : PetscFunctionReturn(PETSC_SUCCESS);
479 : }
480 :
481 16 : static PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
482 : {
483 16 : PetscFunctionBegin;
484 16 : PetscCall(PetscFree(pep->data));
485 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL));
486 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL));
487 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL));
488 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL));
489 16 : PetscFunctionReturn(PETSC_SUCCESS);
490 : }
491 :
492 16 : SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
493 : {
494 16 : PEP_QARNOLDI *ctx;
495 :
496 16 : PetscFunctionBegin;
497 16 : PetscCall(PetscNew(&ctx));
498 16 : pep->data = (void*)ctx;
499 :
500 16 : pep->lineariz = PETSC_TRUE;
501 16 : ctx->lock = PETSC_TRUE;
502 :
503 16 : pep->ops->solve = PEPSolve_QArnoldi;
504 16 : pep->ops->setup = PEPSetUp_QArnoldi;
505 16 : pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
506 16 : pep->ops->destroy = PEPDestroy_QArnoldi;
507 16 : pep->ops->view = PEPView_QArnoldi;
508 16 : pep->ops->backtransform = PEPBackTransform_Default;
509 16 : pep->ops->computevectors = PEPComputeVectors_Default;
510 16 : pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
511 16 : pep->ops->setdefaultst = PEPSetDefaultST_Transform;
512 :
513 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi));
514 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi));
515 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi));
516 16 : PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi));
517 16 : PetscFunctionReturn(PETSC_SUCCESS);
518 : }
|
