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