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 29 : static PetscErrorCode NEPSetUp_SLP(NEP nep)
35 : {
36 29 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
37 29 : PetscBool flg;
38 29 : ST st;
39 :
40 29 : PetscFunctionBegin;
41 29 : if (nep->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
42 29 : nep->ncv = nep->nev;
43 29 : if (nep->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
44 29 : nep->mpd = nep->nev;
45 29 : 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 29 : if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
47 29 : if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
48 29 : PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
49 29 : NEPCheckUnsupported(nep,NEP_FEATURE_REGION);
50 :
51 29 : if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
52 29 : PetscCall(EPSGetST(ctx->eps,&st));
53 29 : PetscCall(PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,""));
54 29 : PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"SLP does not support spectral transformation");
55 29 : PetscCall(EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE));
56 29 : PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_LARGEST_MAGNITUDE));
57 43 : PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(nep->tol)/10.0,nep->max_it));
58 29 : if (nep->tol==(PetscReal)PETSC_DETERMINE) nep->tol = SLEPC_DEFAULT_TOL;
59 29 : if (ctx->deftol==(PetscReal)PETSC_DETERMINE) ctx->deftol = nep->tol;
60 :
61 29 : 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 24 : nep->ops->solve = NEPSolve_SLP;
73 24 : nep->ops->computevectors = NEPComputeVectors_Schur;
74 : }
75 29 : PetscCall(NEPAllocateSolution(nep,0));
76 29 : PetscFunctionReturn(PETSC_SUCCESS);
77 : }
78 :
79 1357 : static PetscErrorCode MatMult_SLP(Mat M,Vec x,Vec y)
80 : {
81 1357 : NEP_SLP_MATSHELL *ctx;
82 :
83 1357 : PetscFunctionBegin;
84 1357 : PetscCall(MatShellGetContext(M,&ctx));
85 1357 : PetscCall(MatMult(ctx->extop->MJ,x,ctx->w));
86 1357 : PetscCall(NEPDeflationFunctionSolve(ctx->extop,ctx->w,y));
87 1357 : PetscFunctionReturn(PETSC_SUCCESS);
88 : }
89 :
90 24 : static PetscErrorCode MatDestroy_SLP(Mat M)
91 : {
92 24 : NEP_SLP_MATSHELL *ctx;
93 :
94 24 : PetscFunctionBegin;
95 24 : PetscCall(MatShellGetContext(M,&ctx));
96 24 : PetscCall(VecDestroy(&ctx->w));
97 24 : PetscCall(PetscFree(ctx));
98 24 : PetscFunctionReturn(PETSC_SUCCESS);
99 : }
100 :
101 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
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 80 : static PetscErrorCode NEPSLPSetUpLinearEP(NEP nep,NEP_EXT_OP extop,PetscScalar lambda,Vec u,PetscBool ini)
115 : {
116 80 : NEP_SLP *slpctx = (NEP_SLP*)nep->data;
117 80 : Mat Mshell;
118 80 : PetscInt nloc,mloc;
119 80 : NEP_SLP_MATSHELL *shellctx;
120 :
121 80 : PetscFunctionBegin;
122 80 : if (ini) {
123 : /* Create mat shell */
124 24 : PetscCall(PetscNew(&shellctx));
125 24 : shellctx->extop = extop;
126 24 : PetscCall(NEPDeflationCreateVec(extop,&shellctx->w));
127 24 : PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
128 24 : nloc += extop->szd; mloc += extop->szd;
129 24 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell));
130 24 : PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLP));
131 24 : PetscCall(MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_SLP));
132 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
133 : PetscCall(MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_SLP));
134 : #endif
135 24 : PetscCall(EPSSetOperators(slpctx->eps,Mshell,NULL));
136 24 : PetscCall(MatDestroy(&Mshell));
137 : }
138 80 : PetscCall(NEPDeflationSolveSetUp(extop,lambda));
139 80 : PetscCall(NEPDeflationComputeJacobian(extop,lambda,NULL));
140 80 : PetscCall(EPSSetInitialSpace(slpctx->eps,1,&u));
141 80 : PetscFunctionReturn(PETSC_SUCCESS);
142 : }
143 :
144 24 : PetscErrorCode NEPSolve_SLP(NEP nep)
145 : {
146 24 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
147 24 : Mat F,H,A;
148 24 : Vec uu,u,r;
149 24 : PetscScalar sigma,lambda,mu,im;
150 24 : PetscReal resnorm;
151 24 : PetscInt nconv;
152 24 : PetscBool skip=PETSC_FALSE,lock=PETSC_FALSE;
153 24 : NEP_EXT_OP extop=NULL; /* Extended operator for deflation */
154 :
155 24 : PetscFunctionBegin;
156 : /* get initial approximation of eigenvalue and eigenvector */
157 24 : PetscCall(NEPGetDefaultShift(nep,&sigma));
158 24 : if (!nep->nini) PetscCall(BVSetRandomColumn(nep->V,0));
159 24 : lambda = sigma;
160 24 : if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
161 24 : PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_TRUE,nep->nev,&extop));
162 24 : PetscCall(NEPDeflationCreateVec(extop,&u));
163 24 : PetscCall(VecDuplicate(u,&r));
164 24 : PetscCall(BVGetColumn(nep->V,0,&uu));
165 24 : PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE));
166 24 : PetscCall(BVRestoreColumn(nep->V,0,&uu));
167 :
168 : /* Restart loop */
169 144 : while (nep->reason == NEP_CONVERGED_ITERATING) {
170 120 : nep->its++;
171 :
172 : /* form residual, r = T(lambda)*u (used in convergence test only) */
173 120 : PetscCall(NEPDeflationComputeFunction(extop,lambda,&F));
174 120 : PetscCall(MatMult(F,u,r));
175 :
176 : /* convergence test */
177 120 : PetscCall(VecNorm(r,NORM_2,&resnorm));
178 120 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
179 120 : nep->eigr[nep->nconv] = lambda;
180 120 : if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
181 40 : if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
182 16 : PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
183 16 : PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
184 16 : PetscCall(MatMult(F,u,r));
185 16 : PetscCall(VecNorm(r,NORM_2,&resnorm));
186 16 : PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
187 16 : if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
188 24 : } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
189 : }
190 :
191 80 : if (lock) {
192 40 : PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
193 40 : nep->nconv = nep->nconv + 1;
194 40 : skip = PETSC_TRUE;
195 40 : lock = PETSC_FALSE;
196 40 : PetscCall(NEPDeflationLocking(extop,u,lambda));
197 : }
198 120 : PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
199 120 : 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 120 : if (nep->reason == NEP_CONVERGED_ITERATING) {
202 96 : if (!skip) {
203 : /* evaluate T(lambda) and T'(lambda) */
204 80 : PetscCall(NEPSLPSetUpLinearEP(nep,extop,lambda,u,nep->its==1?PETSC_TRUE:PETSC_FALSE));
205 : /* compute new eigenvalue correction mu and eigenvector approximation u */
206 80 : PetscCall(EPSSolve(ctx->eps));
207 80 : PetscCall(EPSGetConverged(ctx->eps,&nconv));
208 80 : 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 80 : PetscCall(EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL));
214 80 : mu = 1.0/mu;
215 80 : PetscCheck(PetscAbsScalar(im)<PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
216 : } else {
217 16 : nep->its--; /* do not count this as a full iteration */
218 : /* use second eigenpair computed in previous iteration */
219 16 : PetscCall(EPSGetConverged(ctx->eps,&nconv));
220 16 : if (nconv>=2) {
221 16 : PetscCall(EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL));
222 16 : 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 96 : lambda = lambda - mu;
231 : }
232 : }
233 24 : PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
234 24 : PetscCall(DSSetType(nep->ds,DSNHEP));
235 24 : PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
236 24 : PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
237 24 : PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
238 24 : PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
239 24 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
240 24 : PetscCall(MatDestroy(&H));
241 24 : PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
242 24 : PetscCall(NEPDeflationReset(extop));
243 24 : PetscCall(VecDestroy(&u));
244 24 : PetscCall(VecDestroy(&r));
245 24 : PetscFunctionReturn(PETSC_SUCCESS);
246 : }
247 :
248 25 : static PetscErrorCode NEPSetFromOptions_SLP(NEP nep,PetscOptionItems *PetscOptionsObject)
249 : {
250 25 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
251 25 : PetscBool flg;
252 25 : PetscReal r;
253 :
254 25 : PetscFunctionBegin;
255 25 : PetscOptionsHeadBegin(PetscOptionsObject,"NEP SLP Options");
256 :
257 25 : r = 0.0;
258 25 : PetscCall(PetscOptionsReal("-nep_slp_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPSLPSetDeflationThreshold",ctx->deftol,&r,&flg));
259 25 : if (flg) PetscCall(NEPSLPSetDeflationThreshold(nep,r));
260 :
261 25 : PetscOptionsHeadEnd();
262 :
263 25 : if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
264 25 : PetscCall(EPSSetFromOptions(ctx->eps));
265 25 : if (nep->twosided) {
266 5 : if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
267 5 : PetscCall(EPSSetFromOptions(ctx->epsts));
268 : }
269 25 : if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
270 25 : PetscCall(KSPSetFromOptions(ctx->ksp));
271 25 : 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 || deftol == (PetscReal)PETSC_DECIDE) {
280 0 : ctx->deftol = PETSC_DETERMINE;
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 9 : 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 23 : static PetscErrorCode NEPSLPGetEPS_SLP(NEP nep,EPS *eps)
393 : {
394 23 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
395 :
396 23 : PetscFunctionBegin;
397 23 : if (!ctx->eps) {
398 23 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps));
399 23 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1));
400 23 : PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix));
401 23 : PetscCall(EPSAppendOptionsPrefix(ctx->eps,"nep_slp_"));
402 23 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options));
403 : }
404 23 : *eps = ctx->eps;
405 23 : 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 23 : PetscErrorCode NEPSLPGetEPS(NEP nep,EPS *eps)
425 : {
426 23 : PetscFunctionBegin;
427 23 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
428 23 : PetscAssertPointer(eps,2);
429 23 : PetscUseMethod(nep,"NEPSLPGetEPS_C",(NEP,EPS*),(nep,eps));
430 23 : 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 23 : static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
549 : {
550 23 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
551 :
552 23 : PetscFunctionBegin;
553 23 : if (!ctx->ksp) {
554 23 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
555 23 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
556 23 : PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
557 23 : PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_"));
558 23 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
559 23 : PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
560 35 : PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
561 : }
562 23 : *ksp = ctx->ksp;
563 23 : 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 23 : PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
583 : {
584 23 : PetscFunctionBegin;
585 23 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
586 23 : PetscAssertPointer(ksp,2);
587 23 : PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
588 23 : 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 29 : static PetscErrorCode NEPReset_SLP(NEP nep)
615 : {
616 29 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
617 :
618 29 : PetscFunctionBegin;
619 29 : PetscCall(EPSReset(ctx->eps));
620 29 : if (nep->twosided) PetscCall(EPSReset(ctx->epsts));
621 29 : PetscCall(KSPReset(ctx->ksp));
622 29 : PetscFunctionReturn(PETSC_SUCCESS);
623 : }
624 :
625 25 : static PetscErrorCode NEPDestroy_SLP(NEP nep)
626 : {
627 25 : NEP_SLP *ctx = (NEP_SLP*)nep->data;
628 :
629 25 : PetscFunctionBegin;
630 25 : PetscCall(KSPDestroy(&ctx->ksp));
631 25 : PetscCall(EPSDestroy(&ctx->eps));
632 25 : PetscCall(EPSDestroy(&ctx->epsts));
633 25 : PetscCall(PetscFree(nep->data));
634 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NULL));
635 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NULL));
636 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL));
637 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL));
638 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NULL));
639 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NULL));
640 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL));
641 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL));
642 25 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 25 : SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
646 : {
647 25 : NEP_SLP *ctx;
648 :
649 25 : PetscFunctionBegin;
650 25 : PetscCall(PetscNew(&ctx));
651 25 : nep->data = (void*)ctx;
652 :
653 25 : nep->useds = PETSC_TRUE;
654 25 : ctx->deftol = PETSC_DETERMINE;
655 :
656 25 : nep->ops->solve = NEPSolve_SLP;
657 25 : nep->ops->setup = NEPSetUp_SLP;
658 25 : nep->ops->setfromoptions = NEPSetFromOptions_SLP;
659 25 : nep->ops->reset = NEPReset_SLP;
660 25 : nep->ops->destroy = NEPDestroy_SLP;
661 25 : nep->ops->view = NEPView_SLP;
662 25 : nep->ops->computevectors = NEPComputeVectors_Schur;
663 :
664 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NEPSLPSetDeflationThreshold_SLP));
665 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NEPSLPGetDeflationThreshold_SLP));
666 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP));
667 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP));
668 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NEPSLPSetEPSLeft_SLP));
669 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NEPSLPGetEPSLeft_SLP));
670 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP));
671 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP));
672 25 : PetscFunctionReturn(PETSC_SUCCESS);
673 : }
|
