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 nonlinear eigensolver: "rii"
12 :
13 : Method: Residual inverse iteration
14 :
15 : Algorithm:
16 :
17 : Simple residual inverse iteration with varying shift.
18 :
19 : References:
20 :
21 : [1] A. Neumaier, "Residual inverse iteration for the nonlinear
22 : eigenvalue problem", SIAM J. Numer. Anal. 22(5):914-923, 1985.
23 : */
24 :
25 : #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
26 : #include <../src/nep/impls/nepdefl.h>
27 :
28 : typedef struct {
29 : PetscInt max_inner_it; /* maximum number of Newton iterations */
30 : PetscInt lag; /* interval to rebuild preconditioner */
31 : PetscBool cctol; /* constant correction tolerance */
32 : PetscBool herm; /* whether the Hermitian version of the scalar equation must be used */
33 : PetscReal deftol; /* tolerance for the deflation (threshold) */
34 : KSP ksp; /* linear solver object */
35 : } NEP_RII;
36 :
37 24 : static PetscErrorCode NEPSetUp_RII(NEP nep)
38 : {
39 24 : PetscFunctionBegin;
40 24 : if (nep->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
41 24 : nep->ncv = nep->nev;
42 24 : if (nep->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
43 24 : nep->mpd = nep->nev;
44 24 : if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
45 24 : if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
46 24 : PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
47 24 : NEPCheckUnsupported(nep,NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
48 24 : PetscCall(NEPAllocateSolution(nep,0));
49 24 : PetscCall(NEPSetWorkVecs(nep,2));
50 24 : PetscFunctionReturn(PETSC_SUCCESS);
51 : }
52 :
53 24 : static PetscErrorCode NEPSolve_RII(NEP nep)
54 : {
55 24 : NEP_RII *ctx = (NEP_RII*)nep->data;
56 24 : Mat T,Tp,H,A;
57 24 : Vec uu,u,r,delta,t;
58 24 : PetscScalar lambda,lambda2,sigma,a1,a2,corr;
59 24 : PetscReal nrm,resnorm=1.0,ktol=0.1,perr,rtol;
60 24 : PetscBool skip=PETSC_FALSE,lock=PETSC_FALSE;
61 24 : PetscInt inner_its,its=0;
62 24 : NEP_EXT_OP extop=NULL;
63 24 : KSPConvergedReason kspreason;
64 :
65 24 : PetscFunctionBegin;
66 : /* get initial approximation of eigenvalue and eigenvector */
67 24 : PetscCall(NEPGetDefaultShift(nep,&sigma));
68 24 : lambda = sigma;
69 24 : if (!nep->nini) {
70 22 : PetscCall(BVSetRandomColumn(nep->V,0));
71 22 : PetscCall(BVNormColumn(nep->V,0,NORM_2,&nrm));
72 22 : PetscCall(BVScaleColumn(nep->V,0,1.0/nrm));
73 : }
74 24 : if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
75 24 : PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop));
76 24 : PetscCall(NEPDeflationCreateVec(extop,&u));
77 24 : PetscCall(VecDuplicate(u,&r));
78 24 : PetscCall(VecDuplicate(u,&delta));
79 24 : PetscCall(BVGetColumn(nep->V,0,&uu));
80 24 : PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE));
81 24 : PetscCall(BVRestoreColumn(nep->V,0,&uu));
82 :
83 : /* prepare linear solver */
84 24 : PetscCall(NEPDeflationSolveSetUp(extop,sigma));
85 24 : PetscCall(KSPGetTolerances(ctx->ksp,&rtol,NULL,NULL,NULL));
86 :
87 24 : PetscCall(VecCopy(u,r));
88 24 : PetscCall(NEPDeflationFunctionSolve(extop,r,u));
89 24 : PetscCall(VecNorm(u,NORM_2,&nrm));
90 24 : PetscCall(VecScale(u,1.0/nrm));
91 :
92 : /* Restart loop */
93 381 : while (nep->reason == NEP_CONVERGED_ITERATING) {
94 357 : its++;
95 :
96 : /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
97 : estimate as starting value. */
98 357 : inner_its=0;
99 357 : lambda2 = lambda;
100 827 : do {
101 827 : PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
102 827 : PetscCall(MatMult(T,u,r));
103 827 : if (!ctx->herm) {
104 568 : PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
105 568 : PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
106 568 : if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
107 568 : t = delta;
108 259 : } else t = r;
109 827 : PetscCall(VecDot(t,u,&a1));
110 827 : PetscCall(NEPDeflationComputeJacobian(extop,lambda,&Tp));
111 827 : PetscCall(MatMult(Tp,u,r));
112 827 : if (!ctx->herm) {
113 568 : PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
114 568 : PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
115 568 : if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
116 568 : t = delta;
117 259 : } else t = r;
118 827 : PetscCall(VecDot(t,u,&a2));
119 827 : corr = a1/a2;
120 827 : lambda = lambda - corr;
121 827 : inner_its++;
122 827 : } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);
123 :
124 : /* form residual, r = T(lambda)*u */
125 357 : PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
126 357 : PetscCall(MatMult(T,u,r));
127 :
128 : /* convergence test */
129 357 : perr = nep->errest[nep->nconv];
130 357 : PetscCall(VecNorm(r,NORM_2,&resnorm));
131 357 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
132 357 : nep->eigr[nep->nconv] = lambda;
133 357 : if (its>1 && (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol)) {
134 36 : if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
135 4 : PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
136 4 : PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
137 4 : PetscCall(MatMult(T,u,r));
138 4 : PetscCall(VecNorm(r,NORM_2,&resnorm));
139 4 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
140 4 : if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
141 32 : } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
142 : }
143 322 : if (lock) {
144 35 : PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
145 35 : nep->nconv = nep->nconv + 1;
146 35 : PetscCall(NEPDeflationLocking(extop,u,lambda));
147 35 : nep->its += its;
148 35 : skip = PETSC_TRUE;
149 35 : lock = PETSC_FALSE;
150 : }
151 357 : PetscCall((*nep->stopping)(nep,nep->its+its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
152 357 : if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its+its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));
153 :
154 357 : if (nep->reason == NEP_CONVERGED_ITERATING) {
155 333 : if (!skip) {
156 : /* update preconditioner and set adaptive tolerance */
157 322 : if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) PetscCall(NEPDeflationSolveSetUp(extop,lambda2));
158 322 : if (!ctx->cctol) {
159 314 : ktol = PetscMax(ktol/2.0,rtol);
160 314 : PetscCall(KSPSetTolerances(ctx->ksp,ktol,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
161 : }
162 :
163 : /* eigenvector correction: delta = T(sigma)\r */
164 322 : PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
165 322 : PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
166 322 : if (kspreason<0) {
167 0 : PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its));
168 0 : nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
169 0 : break;
170 : }
171 :
172 : /* update eigenvector: u = u - delta */
173 322 : PetscCall(VecAXPY(u,-1.0,delta));
174 :
175 : /* normalize eigenvector */
176 322 : PetscCall(VecNormalize(u,NULL));
177 : } else {
178 11 : its = -1;
179 11 : PetscCall(NEPDeflationSetRandomVec(extop,u));
180 11 : PetscCall(NEPDeflationSolveSetUp(extop,sigma));
181 11 : PetscCall(VecCopy(u,r));
182 11 : PetscCall(NEPDeflationFunctionSolve(extop,r,u));
183 11 : PetscCall(VecNorm(u,NORM_2,&nrm));
184 11 : PetscCall(VecScale(u,1.0/nrm));
185 11 : lambda = sigma;
186 11 : skip = PETSC_FALSE;
187 : }
188 : }
189 : }
190 24 : PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
191 24 : PetscCall(DSSetType(nep->ds,DSNHEP));
192 24 : PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
193 24 : PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
194 24 : PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
195 24 : PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
196 24 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
197 24 : PetscCall(MatDestroy(&H));
198 24 : PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
199 24 : PetscCall(NEPDeflationReset(extop));
200 24 : PetscCall(VecDestroy(&u));
201 24 : PetscCall(VecDestroy(&r));
202 24 : PetscCall(VecDestroy(&delta));
203 24 : PetscFunctionReturn(PETSC_SUCCESS);
204 : }
205 :
206 20 : static PetscErrorCode NEPSetFromOptions_RII(NEP nep,PetscOptionItems *PetscOptionsObject)
207 : {
208 20 : NEP_RII *ctx = (NEP_RII*)nep->data;
209 20 : PetscBool flg;
210 20 : PetscInt i;
211 20 : PetscReal r;
212 :
213 20 : PetscFunctionBegin;
214 20 : PetscOptionsHeadBegin(PetscOptionsObject,"NEP RII Options");
215 :
216 20 : i = 0;
217 20 : PetscCall(PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&i,&flg));
218 20 : if (flg) PetscCall(NEPRIISetMaximumIterations(nep,i));
219 :
220 20 : PetscCall(PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL));
221 :
222 20 : PetscCall(PetscOptionsBool("-nep_rii_hermitian","Use Hermitian version of the scalar nonlinear equation","NEPRIISetHermitian",ctx->herm,&ctx->herm,NULL));
223 :
224 20 : i = 0;
225 20 : PetscCall(PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg));
226 20 : if (flg) PetscCall(NEPRIISetLagPreconditioner(nep,i));
227 :
228 20 : r = 0.0;
229 20 : PetscCall(PetscOptionsReal("-nep_rii_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPRIISetDeflationThreshold",ctx->deftol,&r,&flg));
230 20 : if (flg) PetscCall(NEPRIISetDeflationThreshold(nep,r));
231 :
232 20 : PetscOptionsHeadEnd();
233 :
234 20 : if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
235 20 : PetscCall(KSPSetFromOptions(ctx->ksp));
236 20 : PetscFunctionReturn(PETSC_SUCCESS);
237 : }
238 :
239 1 : static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
240 : {
241 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
242 :
243 1 : PetscFunctionBegin;
244 1 : if (its==PETSC_DEFAULT || its==PETSC_DECIDE) ctx->max_inner_it = 10;
245 : else {
246 1 : PetscCheck(its>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations must be >0");
247 1 : ctx->max_inner_it = its;
248 : }
249 1 : PetscFunctionReturn(PETSC_SUCCESS);
250 : }
251 :
252 : /*@
253 : NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
254 : used in the RII solver. These are the Newton iterations related to the computation
255 : of the nonlinear Rayleigh functional.
256 :
257 : Logically Collective
258 :
259 : Input Parameters:
260 : + nep - nonlinear eigenvalue solver
261 : - its - maximum inner iterations
262 :
263 : Level: advanced
264 :
265 : .seealso: NEPRIIGetMaximumIterations()
266 : @*/
267 1 : PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
268 : {
269 1 : PetscFunctionBegin;
270 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
271 3 : PetscValidLogicalCollectiveInt(nep,its,2);
272 1 : PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
273 1 : PetscFunctionReturn(PETSC_SUCCESS);
274 : }
275 :
276 1 : static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
277 : {
278 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
279 :
280 1 : PetscFunctionBegin;
281 1 : *its = ctx->max_inner_it;
282 1 : PetscFunctionReturn(PETSC_SUCCESS);
283 : }
284 :
285 : /*@
286 : NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.
287 :
288 : Not Collective
289 :
290 : Input Parameter:
291 : . nep - nonlinear eigenvalue solver
292 :
293 : Output Parameter:
294 : . its - maximum inner iterations
295 :
296 : Level: advanced
297 :
298 : .seealso: NEPRIISetMaximumIterations()
299 : @*/
300 1 : PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
301 : {
302 1 : PetscFunctionBegin;
303 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
304 1 : PetscAssertPointer(its,2);
305 1 : PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
306 1 : PetscFunctionReturn(PETSC_SUCCESS);
307 : }
308 :
309 7 : static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
310 : {
311 7 : NEP_RII *ctx = (NEP_RII*)nep->data;
312 :
313 7 : PetscFunctionBegin;
314 7 : PetscCheck(lag>=0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
315 7 : ctx->lag = lag;
316 7 : PetscFunctionReturn(PETSC_SUCCESS);
317 : }
318 :
319 : /*@
320 : NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
321 : nonlinear solve.
322 :
323 : Logically Collective
324 :
325 : Input Parameters:
326 : + nep - nonlinear eigenvalue solver
327 : - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
328 : computed within the nonlinear iteration, 2 means every second time
329 : the Jacobian is built, etc.
330 :
331 : Options Database Keys:
332 : . -nep_rii_lag_preconditioner <lag> - the lag value
333 :
334 : Notes:
335 : The default is 1.
336 : The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
337 :
338 : Level: intermediate
339 :
340 : .seealso: NEPRIIGetLagPreconditioner()
341 : @*/
342 21 : PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
343 : {
344 21 : PetscFunctionBegin;
345 21 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
346 63 : PetscValidLogicalCollectiveInt(nep,lag,2);
347 21 : PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
348 21 : PetscFunctionReturn(PETSC_SUCCESS);
349 : }
350 :
351 1 : static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
352 : {
353 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
354 :
355 1 : PetscFunctionBegin;
356 1 : *lag = ctx->lag;
357 1 : PetscFunctionReturn(PETSC_SUCCESS);
358 : }
359 :
360 : /*@
361 : NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
362 :
363 : Not Collective
364 :
365 : Input Parameter:
366 : . nep - nonlinear eigenvalue solver
367 :
368 : Output Parameter:
369 : . lag - the lag parameter
370 :
371 : Level: intermediate
372 :
373 : .seealso: NEPRIISetLagPreconditioner()
374 : @*/
375 1 : PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
376 : {
377 1 : PetscFunctionBegin;
378 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
379 1 : PetscAssertPointer(lag,2);
380 1 : PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
381 1 : PetscFunctionReturn(PETSC_SUCCESS);
382 : }
383 :
384 1 : static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
385 : {
386 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
387 :
388 1 : PetscFunctionBegin;
389 1 : ctx->cctol = cct;
390 1 : PetscFunctionReturn(PETSC_SUCCESS);
391 : }
392 :
393 : /*@
394 : NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
395 : in the linear solver constant.
396 :
397 : Logically Collective
398 :
399 : Input Parameters:
400 : + nep - nonlinear eigenvalue solver
401 : - cct - a boolean value
402 :
403 : Options Database Keys:
404 : . -nep_rii_const_correction_tol <bool> - set the boolean flag
405 :
406 : Notes:
407 : By default, an exponentially decreasing tolerance is set in the KSP used
408 : within the nonlinear iteration, so that each Newton iteration requests
409 : better accuracy than the previous one. The constant correction tolerance
410 : flag stops this behaviour.
411 :
412 : Level: intermediate
413 :
414 : .seealso: NEPRIIGetConstCorrectionTol()
415 : @*/
416 1 : PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
417 : {
418 1 : PetscFunctionBegin;
419 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
420 3 : PetscValidLogicalCollectiveBool(nep,cct,2);
421 1 : PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
422 1 : PetscFunctionReturn(PETSC_SUCCESS);
423 : }
424 :
425 1 : static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
426 : {
427 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
428 :
429 1 : PetscFunctionBegin;
430 1 : *cct = ctx->cctol;
431 1 : PetscFunctionReturn(PETSC_SUCCESS);
432 : }
433 :
434 : /*@
435 : NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.
436 :
437 : Not Collective
438 :
439 : Input Parameter:
440 : . nep - nonlinear eigenvalue solver
441 :
442 : Output Parameter:
443 : . cct - the value of the constant tolerance flag
444 :
445 : Level: intermediate
446 :
447 : .seealso: NEPRIISetConstCorrectionTol()
448 : @*/
449 1 : PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
450 : {
451 1 : PetscFunctionBegin;
452 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
453 1 : PetscAssertPointer(cct,2);
454 1 : PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
455 1 : PetscFunctionReturn(PETSC_SUCCESS);
456 : }
457 :
458 1 : static PetscErrorCode NEPRIISetHermitian_RII(NEP nep,PetscBool herm)
459 : {
460 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
461 :
462 1 : PetscFunctionBegin;
463 1 : ctx->herm = herm;
464 1 : PetscFunctionReturn(PETSC_SUCCESS);
465 : }
466 :
467 : /*@
468 : NEPRIISetHermitian - Sets a flag to indicate if the Hermitian version of the
469 : scalar nonlinear equation must be used by the solver.
470 :
471 : Logically Collective
472 :
473 : Input Parameters:
474 : + nep - nonlinear eigenvalue solver
475 : - herm - a boolean value
476 :
477 : Options Database Keys:
478 : . -nep_rii_hermitian <bool> - set the boolean flag
479 :
480 : Notes:
481 : By default, the scalar nonlinear equation x'*inv(T(sigma))*T(z)*x=0 is solved
482 : at each step of the nonlinear iteration. When this flag is set the simpler
483 : form x'*T(z)*x=0 is used, which is supposed to be valid only for Hermitian
484 : problems.
485 :
486 : Level: intermediate
487 :
488 : .seealso: NEPRIIGetHermitian()
489 : @*/
490 1 : PetscErrorCode NEPRIISetHermitian(NEP nep,PetscBool herm)
491 : {
492 1 : PetscFunctionBegin;
493 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
494 3 : PetscValidLogicalCollectiveBool(nep,herm,2);
495 1 : PetscTryMethod(nep,"NEPRIISetHermitian_C",(NEP,PetscBool),(nep,herm));
496 1 : PetscFunctionReturn(PETSC_SUCCESS);
497 : }
498 :
499 1 : static PetscErrorCode NEPRIIGetHermitian_RII(NEP nep,PetscBool *herm)
500 : {
501 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
502 :
503 1 : PetscFunctionBegin;
504 1 : *herm = ctx->herm;
505 1 : PetscFunctionReturn(PETSC_SUCCESS);
506 : }
507 :
508 : /*@
509 : NEPRIIGetHermitian - Returns the flag about using the Hermitian version of
510 : the scalar nonlinear equation.
511 :
512 : Not Collective
513 :
514 : Input Parameter:
515 : . nep - nonlinear eigenvalue solver
516 :
517 : Output Parameter:
518 : . herm - the value of the hermitian flag
519 :
520 : Level: intermediate
521 :
522 : .seealso: NEPRIISetHermitian()
523 : @*/
524 1 : PetscErrorCode NEPRIIGetHermitian(NEP nep,PetscBool *herm)
525 : {
526 1 : PetscFunctionBegin;
527 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
528 1 : PetscAssertPointer(herm,2);
529 1 : PetscUseMethod(nep,"NEPRIIGetHermitian_C",(NEP,PetscBool*),(nep,herm));
530 1 : PetscFunctionReturn(PETSC_SUCCESS);
531 : }
532 :
533 2 : static PetscErrorCode NEPRIISetDeflationThreshold_RII(NEP nep,PetscReal deftol)
534 : {
535 2 : NEP_RII *ctx = (NEP_RII*)nep->data;
536 :
537 2 : PetscFunctionBegin;
538 2 : ctx->deftol = deftol;
539 2 : PetscFunctionReturn(PETSC_SUCCESS);
540 : }
541 :
542 : /*@
543 : NEPRIISetDeflationThreshold - Sets the threshold value used to switch between
544 : deflated and non-deflated iteration.
545 :
546 : Logically Collective
547 :
548 : Input Parameters:
549 : + nep - nonlinear eigenvalue solver
550 : - deftol - the threshold value
551 :
552 : Options Database Keys:
553 : . -nep_rii_deflation_threshold <deftol> - set the threshold
554 :
555 : Notes:
556 : Normally, the solver iterates on the extended problem in order to deflate
557 : previously converged eigenpairs. If this threshold is set to a nonzero value,
558 : then once the residual error is below this threshold the solver will
559 : continue the iteration without deflation. The intention is to be able to
560 : improve the current eigenpair further, despite having previous eigenpairs
561 : with somewhat bad precision.
562 :
563 : Level: advanced
564 :
565 : .seealso: NEPRIIGetDeflationThreshold()
566 : @*/
567 2 : PetscErrorCode NEPRIISetDeflationThreshold(NEP nep,PetscReal deftol)
568 : {
569 2 : PetscFunctionBegin;
570 2 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
571 6 : PetscValidLogicalCollectiveReal(nep,deftol,2);
572 2 : PetscTryMethod(nep,"NEPRIISetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
573 2 : PetscFunctionReturn(PETSC_SUCCESS);
574 : }
575 :
576 1 : static PetscErrorCode NEPRIIGetDeflationThreshold_RII(NEP nep,PetscReal *deftol)
577 : {
578 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
579 :
580 1 : PetscFunctionBegin;
581 1 : *deftol = ctx->deftol;
582 1 : PetscFunctionReturn(PETSC_SUCCESS);
583 : }
584 :
585 : /*@
586 : NEPRIIGetDeflationThreshold - Returns the threshold value that controls deflation.
587 :
588 : Not Collective
589 :
590 : Input Parameter:
591 : . nep - nonlinear eigenvalue solver
592 :
593 : Output Parameter:
594 : . deftol - the threshold
595 :
596 : Level: advanced
597 :
598 : .seealso: NEPRIISetDeflationThreshold()
599 : @*/
600 1 : PetscErrorCode NEPRIIGetDeflationThreshold(NEP nep,PetscReal *deftol)
601 : {
602 1 : PetscFunctionBegin;
603 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
604 1 : PetscAssertPointer(deftol,2);
605 1 : PetscUseMethod(nep,"NEPRIIGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
606 1 : PetscFunctionReturn(PETSC_SUCCESS);
607 : }
608 :
609 1 : static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
610 : {
611 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
612 :
613 1 : PetscFunctionBegin;
614 1 : PetscCall(PetscObjectReference((PetscObject)ksp));
615 1 : PetscCall(KSPDestroy(&ctx->ksp));
616 1 : ctx->ksp = ksp;
617 1 : nep->state = NEP_STATE_INITIAL;
618 1 : PetscFunctionReturn(PETSC_SUCCESS);
619 : }
620 :
621 : /*@
622 : NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
623 : eigenvalue solver.
624 :
625 : Collective
626 :
627 : Input Parameters:
628 : + nep - eigenvalue solver
629 : - ksp - the linear solver object
630 :
631 : Level: advanced
632 :
633 : .seealso: NEPRIIGetKSP()
634 : @*/
635 1 : PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
636 : {
637 1 : PetscFunctionBegin;
638 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
639 1 : PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
640 1 : PetscCheckSameComm(nep,1,ksp,2);
641 1 : PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
642 1 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 22 : static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
646 : {
647 22 : NEP_RII *ctx = (NEP_RII*)nep->data;
648 :
649 22 : PetscFunctionBegin;
650 22 : if (!ctx->ksp) {
651 22 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
652 22 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
653 22 : PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
654 22 : PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_"));
655 22 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
656 22 : PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
657 31 : PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
658 : }
659 22 : *ksp = ctx->ksp;
660 22 : PetscFunctionReturn(PETSC_SUCCESS);
661 : }
662 :
663 : /*@
664 : NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
665 : the nonlinear eigenvalue solver.
666 :
667 : Collective
668 :
669 : Input Parameter:
670 : . nep - nonlinear eigenvalue solver
671 :
672 : Output Parameter:
673 : . ksp - the linear solver object
674 :
675 : Level: advanced
676 :
677 : .seealso: NEPRIISetKSP()
678 : @*/
679 22 : PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
680 : {
681 22 : PetscFunctionBegin;
682 22 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
683 22 : PetscAssertPointer(ksp,2);
684 22 : PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
685 22 : PetscFunctionReturn(PETSC_SUCCESS);
686 : }
687 :
688 1 : static PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
689 : {
690 1 : NEP_RII *ctx = (NEP_RII*)nep->data;
691 1 : PetscBool isascii;
692 :
693 1 : PetscFunctionBegin;
694 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
695 1 : if (isascii) {
696 1 : PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of inner iterations: %" PetscInt_FMT "\n",ctx->max_inner_it));
697 1 : if (ctx->cctol) PetscCall(PetscViewerASCIIPrintf(viewer," using a constant tolerance for the linear solver\n"));
698 1 : if (ctx->herm) PetscCall(PetscViewerASCIIPrintf(viewer," using the Hermitian version of the scalar nonlinear equation\n"));
699 1 : if (ctx->lag) PetscCall(PetscViewerASCIIPrintf(viewer," updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag));
700 1 : if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer," deflation threshold: %g\n",(double)ctx->deftol));
701 1 : if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
702 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
703 1 : PetscCall(KSPView(ctx->ksp,viewer));
704 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
705 : }
706 1 : PetscFunctionReturn(PETSC_SUCCESS);
707 : }
708 :
709 24 : static PetscErrorCode NEPReset_RII(NEP nep)
710 : {
711 24 : NEP_RII *ctx = (NEP_RII*)nep->data;
712 :
713 24 : PetscFunctionBegin;
714 24 : PetscCall(KSPReset(ctx->ksp));
715 24 : PetscFunctionReturn(PETSC_SUCCESS);
716 : }
717 :
718 23 : static PetscErrorCode NEPDestroy_RII(NEP nep)
719 : {
720 23 : NEP_RII *ctx = (NEP_RII*)nep->data;
721 :
722 23 : PetscFunctionBegin;
723 23 : PetscCall(KSPDestroy(&ctx->ksp));
724 23 : PetscCall(PetscFree(nep->data));
725 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL));
726 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL));
727 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL));
728 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL));
729 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL));
730 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL));
731 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NULL));
732 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NULL));
733 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NULL));
734 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NULL));
735 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL));
736 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL));
737 23 : PetscFunctionReturn(PETSC_SUCCESS);
738 : }
739 :
740 23 : SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
741 : {
742 23 : NEP_RII *ctx;
743 :
744 23 : PetscFunctionBegin;
745 23 : PetscCall(PetscNew(&ctx));
746 23 : nep->data = (void*)ctx;
747 23 : ctx->max_inner_it = 10;
748 23 : ctx->lag = 1;
749 23 : ctx->cctol = PETSC_FALSE;
750 23 : ctx->herm = PETSC_FALSE;
751 23 : ctx->deftol = 0.0;
752 :
753 23 : nep->useds = PETSC_TRUE;
754 :
755 23 : nep->ops->solve = NEPSolve_RII;
756 23 : nep->ops->setup = NEPSetUp_RII;
757 23 : nep->ops->setfromoptions = NEPSetFromOptions_RII;
758 23 : nep->ops->reset = NEPReset_RII;
759 23 : nep->ops->destroy = NEPDestroy_RII;
760 23 : nep->ops->view = NEPView_RII;
761 23 : nep->ops->computevectors = NEPComputeVectors_Schur;
762 :
763 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII));
764 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII));
765 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII));
766 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII));
767 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII));
768 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII));
769 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NEPRIISetHermitian_RII));
770 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NEPRIIGetHermitian_RII));
771 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NEPRIISetDeflationThreshold_RII));
772 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NEPRIIGetDeflationThreshold_RII));
773 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII));
774 23 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII));
775 23 : PetscFunctionReturn(PETSC_SUCCESS);
776 : }
|
