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: "slp"
12 :
13 : Method: Successive linear problems
14 :
15 : Algorithm:
16 :
17 : Newton-type iteration based on first order Taylor approximation.
18 :
19 : References:
20 :
21 : [1] A. Ruhe, "Algorithms for the nonlinear eigenvalue problem", SIAM J.
22 : Numer. Anal. 10(4):674-689, 1973.
23 : */
24 :
25 : #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
26 : #include <../src/nep/impls/nepdefl.h>
27 : #include "slp.h"
28 :
29 : typedef struct {
30 : NEP_EXT_OP extop;
31 : Vec w;
32 : } NEP_SLP_MATSHELL;
33 :
34 30 : static PetscErrorCode NEPSetUp_SLP(NEP nep)
35 : {
36 30 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
37 30 : PetscBool flg;
38 30 : ST st;
39 :
40 30 : PetscFunctionBegin;
41 30 : if (nep->ncv!=PETSC_DEFAULT) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
42 30 : nep->ncv = nep->nev;
43 30 : if (nep->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
44 30 : nep->mpd = nep->nev;
45 30 : PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
46 30 : if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
47 30 : if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
48 30 : PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
49 30 : NEPCheckUnsupported(nep,NEP_FEATURE_REGION);
50 :
51 30 : if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
52 30 : PetscCall(EPSGetST(ctx->eps,&st));
53 30 : PetscCall(PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,""));
54 30 : PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"SLP does not support spectral transformation");
55 30 : PetscCall(EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE));
56 30 : PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_LARGEST_MAGNITUDE));
57 45 : PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(nep->tol)/10.0,nep->max_it));
58 30 : if (nep->tol==(PetscReal)PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
59 30 : if (ctx->deftol==(PetscReal)PETSC_DEFAULT) ctx->deftol = nep->tol;
60 :
61 30 : if (nep->twosided) {
62 5 : nep->ops->solve = NEPSolve_SLP_Twosided;
63 5 : nep->ops->computevectors = NULL;
64 5 : if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
65 5 : PetscCall(EPSGetST(ctx->epsts,&st));
66 5 : PetscCall(PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,""));
67 5 : PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"SLP does not support spectral transformation");
68 5 : PetscCall(EPSSetDimensions(ctx->epsts,1,PETSC_DECIDE,PETSC_DECIDE));
69 5 : PetscCall(EPSSetWhichEigenpairs(ctx->epsts,EPS_LARGEST_MAGNITUDE));
70 5 : PetscCall(EPSSetTolerances(ctx->epsts,SlepcDefaultTol(nep->tol)/10.0,nep->max_it));
71 : } else {
72 25 : nep->ops->solve = NEPSolve_SLP;
73 25 : nep->ops->computevectors = NEPComputeVectors_Schur;
74 : }
75 30 : PetscCall(NEPAllocateSolution(nep,0));
76 30 : PetscFunctionReturn(PETSC_SUCCESS);
77 : }
78 :
79 1454 : static PetscErrorCode MatMult_SLP(Mat M,Vec x,Vec y)
80 : {
81 1454 : NEP_SLP_MATSHELL *ctx;
82 :
83 1454 : PetscFunctionBegin;
84 1454 : PetscCall(MatShellGetContext(M,&ctx));
85 1454 : PetscCall(MatMult(ctx->extop->MJ,x,ctx->w));
86 1454 : PetscCall(NEPDeflationFunctionSolve(ctx->extop,ctx->w,y));
87 1454 : PetscFunctionReturn(PETSC_SUCCESS);
88 : }
89 :
90 25 : static PetscErrorCode MatDestroy_SLP(Mat M)
91 : {
92 25 : NEP_SLP_MATSHELL *ctx;
93 :
94 25 : PetscFunctionBegin;
95 25 : PetscCall(MatShellGetContext(M,&ctx));
96 25 : PetscCall(VecDestroy(&ctx->w));
97 25 : PetscCall(PetscFree(ctx));
98 25 : PetscFunctionReturn(PETSC_SUCCESS);
99 : }
100 :
101 : #if defined(PETSC_HAVE_CUDA)
102 : static PetscErrorCode MatCreateVecs_SLP(Mat M,Vec *left,Vec *right)
103 : {
104 : NEP_SLP_MATSHELL *ctx;
105 :
106 : PetscFunctionBegin;
107 : PetscCall(MatShellGetContext(M,&ctx));
108 : if (right) PetscCall(VecDuplicate(ctx->w,right));
109 : if (left) PetscCall(VecDuplicate(ctx->w,left));
110 : PetscFunctionReturn(PETSC_SUCCESS);
111 : }
112 : #endif
113 :
114 91 : static PetscErrorCode NEPSLPSetUpLinearEP(NEP nep,NEP_EXT_OP extop,PetscScalar lambda,Vec u,PetscBool ini)
115 : {
116 91 : NEP_SLP *slpctx = (NEP_SLP*)nep->data;
117 91 : Mat Mshell;
118 91 : PetscInt nloc,mloc;
119 91 : NEP_SLP_MATSHELL *shellctx;
120 :
121 91 : PetscFunctionBegin;
122 91 : if (ini) {
123 : /* Create mat shell */
124 25 : PetscCall(PetscNew(&shellctx));
125 25 : shellctx->extop = extop;
126 25 : PetscCall(NEPDeflationCreateVec(extop,&shellctx->w));
127 25 : PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
128 25 : nloc += extop->szd; mloc += extop->szd;
129 25 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell));
130 25 : PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLP));
131 25 : PetscCall(MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_SLP));
132 : #if defined(PETSC_HAVE_CUDA)
133 : PetscCall(MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_SLP));
134 : #endif
135 25 : PetscCall(EPSSetOperators(slpctx->eps,Mshell,NULL));
136 25 : PetscCall(MatDestroy(&Mshell));
137 : }
138 91 : PetscCall(NEPDeflationSolveSetUp(extop,lambda));
139 91 : PetscCall(NEPDeflationComputeJacobian(extop,lambda,NULL));
140 91 : PetscCall(EPSSetInitialSpace(slpctx->eps,1,&u));
141 91 : PetscFunctionReturn(PETSC_SUCCESS);
142 : }
143 :
144 25 : PetscErrorCode NEPSolve_SLP(NEP nep)
145 : {
146 25 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
147 25 : Mat F,H,A;
148 25 : Vec uu,u,r;
149 25 : PetscScalar sigma,lambda,mu,im;
150 25 : PetscReal resnorm;
151 25 : PetscInt nconv;
152 25 : PetscBool skip=PETSC_FALSE,lock=PETSC_FALSE;
153 25 : NEP_EXT_OP extop=NULL; /* Extended operator for deflation */
154 :
155 25 : PetscFunctionBegin;
156 : /* get initial approximation of eigenvalue and eigenvector */
157 25 : PetscCall(NEPGetDefaultShift(nep,&sigma));
158 25 : if (!nep->nini) PetscCall(BVSetRandomColumn(nep->V,0));
159 25 : lambda = sigma;
160 25 : if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
161 25 : PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_TRUE,nep->nev,&extop));
162 25 : PetscCall(NEPDeflationCreateVec(extop,&u));
163 25 : PetscCall(VecDuplicate(u,&r));
164 25 : PetscCall(BVGetColumn(nep->V,0,&uu));
165 25 : PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE));
166 25 : PetscCall(BVRestoreColumn(nep->V,0,&uu));
167 :
168 : /* Restart loop */
169 160 : while (nep->reason == NEP_CONVERGED_ITERATING) {
170 135 : nep->its++;
171 :
172 : /* form residual, r = T(lambda)*u (used in convergence test only) */
173 135 : PetscCall(NEPDeflationComputeFunction(extop,lambda,&F));
174 135 : PetscCall(MatMult(F,u,r));
175 :
176 : /* convergence test */
177 135 : PetscCall(VecNorm(r,NORM_2,&resnorm));
178 135 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
179 135 : nep->eigr[nep->nconv] = lambda;
180 135 : if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
181 44 : if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
182 19 : PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
183 19 : PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
184 19 : PetscCall(MatMult(F,u,r));
185 19 : PetscCall(VecNorm(r,NORM_2,&resnorm));
186 19 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
187 19 : if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
188 25 : } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
189 : }
190 :
191 91 : if (lock) {
192 44 : PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
193 44 : nep->nconv = nep->nconv + 1;
194 44 : skip = PETSC_TRUE;
195 44 : lock = PETSC_FALSE;
196 44 : PetscCall(NEPDeflationLocking(extop,u,lambda));
197 : }
198 135 : PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
199 135 : if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));
200 :
201 135 : if (nep->reason == NEP_CONVERGED_ITERATING) {
202 110 : if (!skip) {
203 : /* evaluate T(lambda) and T'(lambda) */
204 91 : PetscCall(NEPSLPSetUpLinearEP(nep,extop,lambda,u,nep->its==1?PETSC_TRUE:PETSC_FALSE));
205 : /* compute new eigenvalue correction mu and eigenvector approximation u */
206 91 : PetscCall(EPSSolve(ctx->eps));
207 91 : PetscCall(EPSGetConverged(ctx->eps,&nconv));
208 91 : if (!nconv) {
209 0 : PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
210 0 : nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
211 0 : break;
212 : }
213 91 : PetscCall(EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL));
214 91 : mu = 1.0/mu;
215 91 : PetscCheck(PetscAbsScalar(im)<PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
216 : } else {
217 19 : nep->its--; /* do not count this as a full iteration */
218 : /* use second eigenpair computed in previous iteration */
219 19 : PetscCall(EPSGetConverged(ctx->eps,&nconv));
220 19 : if (nconv>=2) {
221 19 : PetscCall(EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL));
222 19 : mu = 1.0/mu;
223 : } else {
224 0 : PetscCall(NEPDeflationSetRandomVec(extop,u));
225 0 : mu = lambda-sigma;
226 : }
227 : skip = PETSC_FALSE;
228 : }
229 : /* correct eigenvalue */
230 110 : lambda = lambda - mu;
231 : }
232 : }
233 25 : PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
234 25 : PetscCall(DSSetType(nep->ds,DSNHEP));
235 25 : PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
236 25 : PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
237 25 : PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
238 25 : PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
239 25 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
240 25 : PetscCall(MatDestroy(&H));
241 25 : PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
242 25 : PetscCall(NEPDeflationReset(extop));
243 25 : PetscCall(VecDestroy(&u));
244 25 : PetscCall(VecDestroy(&r));
245 25 : PetscFunctionReturn(PETSC_SUCCESS);
246 : }
247 :
248 26 : static PetscErrorCode NEPSetFromOptions_SLP(NEP nep,PetscOptionItems *PetscOptionsObject)
249 : {
250 26 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
251 26 : PetscBool flg;
252 26 : PetscReal r;
253 :
254 26 : PetscFunctionBegin;
255 26 : PetscOptionsHeadBegin(PetscOptionsObject,"NEP SLP Options");
256 :
257 26 : r = 0.0;
258 26 : PetscCall(PetscOptionsReal("-nep_slp_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPSLPSetDeflationThreshold",ctx->deftol,&r,&flg));
259 26 : if (flg) PetscCall(NEPSLPSetDeflationThreshold(nep,r));
260 :
261 26 : PetscOptionsHeadEnd();
262 :
263 26 : if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
264 26 : PetscCall(EPSSetFromOptions(ctx->eps));
265 26 : if (nep->twosided) {
266 5 : if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
267 5 : PetscCall(EPSSetFromOptions(ctx->epsts));
268 : }
269 26 : if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
270 26 : PetscCall(KSPSetFromOptions(ctx->ksp));
271 26 : PetscFunctionReturn(PETSC_SUCCESS);
272 : }
273 :
274 3 : static PetscErrorCode NEPSLPSetDeflationThreshold_SLP(NEP nep,PetscReal deftol)
275 : {
276 3 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
277 :
278 3 : PetscFunctionBegin;
279 3 : if (deftol == (PetscReal)PETSC_DEFAULT) {
280 0 : ctx->deftol = PETSC_DEFAULT;
281 0 : nep->state = NEP_STATE_INITIAL;
282 : } else {
283 3 : PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
284 3 : ctx->deftol = deftol;
285 : }
286 3 : PetscFunctionReturn(PETSC_SUCCESS);
287 : }
288 :
289 : /*@
290 : NEPSLPSetDeflationThreshold - Sets the threshold value used to switch between
291 : deflated and non-deflated iteration.
292 :
293 : Logically Collective
294 :
295 : Input Parameters:
296 : + nep - nonlinear eigenvalue solver
297 : - deftol - the threshold value
298 :
299 : Options Database Keys:
300 : . -nep_slp_deflation_threshold <deftol> - set the threshold
301 :
302 : Notes:
303 : Normally, the solver iterates on the extended problem in order to deflate
304 : previously converged eigenpairs. If this threshold is set to a nonzero value,
305 : then once the residual error is below this threshold the solver will
306 : continue the iteration without deflation. The intention is to be able to
307 : improve the current eigenpair further, despite having previous eigenpairs
308 : with somewhat bad precision.
309 :
310 : Level: advanced
311 :
312 : .seealso: NEPSLPGetDeflationThreshold()
313 : @*/
314 3 : PetscErrorCode NEPSLPSetDeflationThreshold(NEP nep,PetscReal deftol)
315 : {
316 3 : PetscFunctionBegin;
317 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
318 12 : PetscValidLogicalCollectiveReal(nep,deftol,2);
319 3 : PetscTryMethod(nep,"NEPSLPSetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
320 3 : PetscFunctionReturn(PETSC_SUCCESS);
321 : }
322 :
323 1 : static PetscErrorCode NEPSLPGetDeflationThreshold_SLP(NEP nep,PetscReal *deftol)
324 : {
325 1 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
326 :
327 1 : PetscFunctionBegin;
328 1 : *deftol = ctx->deftol;
329 1 : PetscFunctionReturn(PETSC_SUCCESS);
330 : }
331 :
332 : /*@
333 : NEPSLPGetDeflationThreshold - Returns the threshold value that controls deflation.
334 :
335 : Not Collective
336 :
337 : Input Parameter:
338 : . nep - nonlinear eigenvalue solver
339 :
340 : Output Parameter:
341 : . deftol - the threshold
342 :
343 : Level: advanced
344 :
345 : .seealso: NEPSLPSetDeflationThreshold()
346 : @*/
347 1 : PetscErrorCode NEPSLPGetDeflationThreshold(NEP nep,PetscReal *deftol)
348 : {
349 1 : PetscFunctionBegin;
350 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
351 1 : PetscAssertPointer(deftol,2);
352 1 : PetscUseMethod(nep,"NEPSLPGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
353 1 : PetscFunctionReturn(PETSC_SUCCESS);
354 : }
355 :
356 2 : static PetscErrorCode NEPSLPSetEPS_SLP(NEP nep,EPS eps)
357 : {
358 2 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
359 :
360 2 : PetscFunctionBegin;
361 2 : PetscCall(PetscObjectReference((PetscObject)eps));
362 2 : PetscCall(EPSDestroy(&ctx->eps));
363 2 : ctx->eps = eps;
364 2 : nep->state = NEP_STATE_INITIAL;
365 2 : PetscFunctionReturn(PETSC_SUCCESS);
366 : }
367 :
368 : /*@
369 : NEPSLPSetEPS - Associate a linear eigensolver object (EPS) to the
370 : nonlinear eigenvalue solver.
371 :
372 : Collective
373 :
374 : Input Parameters:
375 : + nep - nonlinear eigenvalue solver
376 : - eps - the eigensolver object
377 :
378 : Level: advanced
379 :
380 : .seealso: NEPSLPGetEPS()
381 : @*/
382 2 : PetscErrorCode NEPSLPSetEPS(NEP nep,EPS eps)
383 : {
384 2 : PetscFunctionBegin;
385 2 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
386 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
387 2 : PetscCheckSameComm(nep,1,eps,2);
388 2 : PetscTryMethod(nep,"NEPSLPSetEPS_C",(NEP,EPS),(nep,eps));
389 2 : PetscFunctionReturn(PETSC_SUCCESS);
390 : }
391 :
392 24 : static PetscErrorCode NEPSLPGetEPS_SLP(NEP nep,EPS *eps)
393 : {
394 24 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
395 :
396 24 : PetscFunctionBegin;
397 24 : if (!ctx->eps) {
398 24 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps));
399 24 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1));
400 24 : PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix));
401 24 : PetscCall(EPSAppendOptionsPrefix(ctx->eps,"nep_slp_"));
402 24 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options));
403 : }
404 24 : *eps = ctx->eps;
405 24 : PetscFunctionReturn(PETSC_SUCCESS);
406 : }
407 :
408 : /*@
409 : NEPSLPGetEPS - Retrieve the linear eigensolver object (EPS) associated
410 : to the nonlinear eigenvalue solver.
411 :
412 : Collective
413 :
414 : Input Parameter:
415 : . nep - nonlinear eigenvalue solver
416 :
417 : Output Parameter:
418 : . eps - the eigensolver object
419 :
420 : Level: advanced
421 :
422 : .seealso: NEPSLPSetEPS()
423 : @*/
424 24 : PetscErrorCode NEPSLPGetEPS(NEP nep,EPS *eps)
425 : {
426 24 : PetscFunctionBegin;
427 24 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
428 24 : PetscAssertPointer(eps,2);
429 24 : PetscUseMethod(nep,"NEPSLPGetEPS_C",(NEP,EPS*),(nep,eps));
430 24 : PetscFunctionReturn(PETSC_SUCCESS);
431 : }
432 :
433 0 : static PetscErrorCode NEPSLPSetEPSLeft_SLP(NEP nep,EPS eps)
434 : {
435 0 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
436 :
437 0 : PetscFunctionBegin;
438 0 : PetscCall(PetscObjectReference((PetscObject)eps));
439 0 : PetscCall(EPSDestroy(&ctx->epsts));
440 0 : ctx->epsts = eps;
441 0 : nep->state = NEP_STATE_INITIAL;
442 0 : PetscFunctionReturn(PETSC_SUCCESS);
443 : }
444 :
445 : /*@
446 : NEPSLPSetEPSLeft - Associate a linear eigensolver object (EPS) to the
447 : nonlinear eigenvalue solver, used to compute left eigenvectors in the
448 : two-sided variant of SLP.
449 :
450 : Collective
451 :
452 : Input Parameters:
453 : + nep - nonlinear eigenvalue solver
454 : - eps - the eigensolver object
455 :
456 : Level: advanced
457 :
458 : .seealso: NEPSLPGetEPSLeft(), NEPSetTwoSided()
459 : @*/
460 0 : PetscErrorCode NEPSLPSetEPSLeft(NEP nep,EPS eps)
461 : {
462 0 : PetscFunctionBegin;
463 0 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
464 0 : PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
465 0 : PetscCheckSameComm(nep,1,eps,2);
466 0 : PetscTryMethod(nep,"NEPSLPSetEPSLeft_C",(NEP,EPS),(nep,eps));
467 0 : PetscFunctionReturn(PETSC_SUCCESS);
468 : }
469 :
470 5 : static PetscErrorCode NEPSLPGetEPSLeft_SLP(NEP nep,EPS *eps)
471 : {
472 5 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
473 :
474 5 : PetscFunctionBegin;
475 5 : if (!ctx->epsts) {
476 5 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->epsts));
477 5 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->epsts,(PetscObject)nep,1));
478 5 : PetscCall(EPSSetOptionsPrefix(ctx->epsts,((PetscObject)nep)->prefix));
479 5 : PetscCall(EPSAppendOptionsPrefix(ctx->epsts,"nep_slp_left_"));
480 5 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->epsts,((PetscObject)nep)->options));
481 : }
482 5 : *eps = ctx->epsts;
483 5 : PetscFunctionReturn(PETSC_SUCCESS);
484 : }
485 :
486 : /*@
487 : NEPSLPGetEPSLeft - Retrieve the linear eigensolver object (EPS) associated
488 : to the nonlinear eigenvalue solver, used to compute left eigenvectors in the
489 : two-sided variant of SLP.
490 :
491 : Collective
492 :
493 : Input Parameter:
494 : . nep - nonlinear eigenvalue solver
495 :
496 : Output Parameter:
497 : . eps - the eigensolver object
498 :
499 : Level: advanced
500 :
501 : .seealso: NEPSLPSetEPSLeft(), NEPSetTwoSided()
502 : @*/
503 5 : PetscErrorCode NEPSLPGetEPSLeft(NEP nep,EPS *eps)
504 : {
505 5 : PetscFunctionBegin;
506 5 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
507 5 : PetscAssertPointer(eps,2);
508 5 : PetscUseMethod(nep,"NEPSLPGetEPSLeft_C",(NEP,EPS*),(nep,eps));
509 5 : PetscFunctionReturn(PETSC_SUCCESS);
510 : }
511 :
512 2 : static PetscErrorCode NEPSLPSetKSP_SLP(NEP nep,KSP ksp)
513 : {
514 2 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
515 :
516 2 : PetscFunctionBegin;
517 2 : PetscCall(PetscObjectReference((PetscObject)ksp));
518 2 : PetscCall(KSPDestroy(&ctx->ksp));
519 2 : ctx->ksp = ksp;
520 2 : nep->state = NEP_STATE_INITIAL;
521 2 : PetscFunctionReturn(PETSC_SUCCESS);
522 : }
523 :
524 : /*@
525 : NEPSLPSetKSP - Associate a linear solver object (KSP) to the nonlinear
526 : eigenvalue solver.
527 :
528 : Collective
529 :
530 : Input Parameters:
531 : + nep - eigenvalue solver
532 : - ksp - the linear solver object
533 :
534 : Level: advanced
535 :
536 : .seealso: NEPSLPGetKSP()
537 : @*/
538 2 : PetscErrorCode NEPSLPSetKSP(NEP nep,KSP ksp)
539 : {
540 2 : PetscFunctionBegin;
541 2 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
542 2 : PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
543 2 : PetscCheckSameComm(nep,1,ksp,2);
544 2 : PetscTryMethod(nep,"NEPSLPSetKSP_C",(NEP,KSP),(nep,ksp));
545 2 : PetscFunctionReturn(PETSC_SUCCESS);
546 : }
547 :
548 24 : static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
549 : {
550 24 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
551 :
552 24 : PetscFunctionBegin;
553 24 : if (!ctx->ksp) {
554 24 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
555 24 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
556 24 : PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
557 24 : PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_"));
558 24 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
559 24 : PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
560 37 : PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
561 : }
562 24 : *ksp = ctx->ksp;
563 24 : PetscFunctionReturn(PETSC_SUCCESS);
564 : }
565 :
566 : /*@
567 : NEPSLPGetKSP - Retrieve the linear solver object (KSP) associated with
568 : the nonlinear eigenvalue solver.
569 :
570 : Collective
571 :
572 : Input Parameter:
573 : . nep - nonlinear eigenvalue solver
574 :
575 : Output Parameter:
576 : . ksp - the linear solver object
577 :
578 : Level: advanced
579 :
580 : .seealso: NEPSLPSetKSP()
581 : @*/
582 24 : PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
583 : {
584 24 : PetscFunctionBegin;
585 24 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
586 24 : PetscAssertPointer(ksp,2);
587 24 : PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
588 24 : PetscFunctionReturn(PETSC_SUCCESS);
589 : }
590 :
591 1 : static PetscErrorCode NEPView_SLP(NEP nep,PetscViewer viewer)
592 : {
593 1 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
594 1 : PetscBool isascii;
595 :
596 1 : PetscFunctionBegin;
597 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
598 1 : if (isascii) {
599 1 : if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer," deflation threshold: %g\n",(double)ctx->deftol));
600 1 : if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
601 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
602 1 : PetscCall(EPSView(ctx->eps,viewer));
603 1 : if (nep->twosided) {
604 0 : if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
605 0 : PetscCall(EPSView(ctx->epsts,viewer));
606 : }
607 1 : if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
608 1 : PetscCall(KSPView(ctx->ksp,viewer));
609 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
610 : }
611 1 : PetscFunctionReturn(PETSC_SUCCESS);
612 : }
613 :
614 30 : static PetscErrorCode NEPReset_SLP(NEP nep)
615 : {
616 30 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
617 :
618 30 : PetscFunctionBegin;
619 30 : PetscCall(EPSReset(ctx->eps));
620 30 : if (nep->twosided) PetscCall(EPSReset(ctx->epsts));
621 30 : PetscCall(KSPReset(ctx->ksp));
622 30 : PetscFunctionReturn(PETSC_SUCCESS);
623 : }
624 :
625 26 : static PetscErrorCode NEPDestroy_SLP(NEP nep)
626 : {
627 26 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
628 :
629 26 : PetscFunctionBegin;
630 26 : PetscCall(KSPDestroy(&ctx->ksp));
631 26 : PetscCall(EPSDestroy(&ctx->eps));
632 26 : PetscCall(EPSDestroy(&ctx->epsts));
633 26 : PetscCall(PetscFree(nep->data));
634 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NULL));
635 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NULL));
636 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL));
637 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL));
638 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NULL));
639 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NULL));
640 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL));
641 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL));
642 26 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 26 : SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
646 : {
647 26 : NEP_SLP *ctx;
648 :
649 26 : PetscFunctionBegin;
650 26 : PetscCall(PetscNew(&ctx));
651 26 : nep->data = (void*)ctx;
652 :
653 26 : nep->useds = PETSC_TRUE;
654 26 : ctx->deftol = PETSC_DEFAULT;
655 :
656 26 : nep->ops->solve = NEPSolve_SLP;
657 26 : nep->ops->setup = NEPSetUp_SLP;
658 26 : nep->ops->setfromoptions = NEPSetFromOptions_SLP;
659 26 : nep->ops->reset = NEPReset_SLP;
660 26 : nep->ops->destroy = NEPDestroy_SLP;
661 26 : nep->ops->view = NEPView_SLP;
662 26 : nep->ops->computevectors = NEPComputeVectors_Schur;
663 :
664 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NEPSLPSetDeflationThreshold_SLP));
665 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NEPSLPGetDeflationThreshold_SLP));
666 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP));
667 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP));
668 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NEPSLPSetEPSLeft_SLP));
669 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NEPSLPGetEPSLeft_SLP));
670 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP));
671 26 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP));
672 26 : PetscFunctionReturn(PETSC_SUCCESS);
673 : }
|
