Actual source code: slp.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
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"
13: Method: Successive linear problems
15: Algorithm:
17: Newton-type iteration based on first order Taylor approximation.
19: References:
21: [1] A. Ruhe, "Algorithms for the nonlinear eigenvalue problem", SIAM J.
22: Numer. Anal. 10(4):674-689, 1973.
23: */
25: #include <slepc/private/nepimpl.h>
26: #include <../src/nep/impls/nepdefl.h>
27: #include "slp.h"
29: typedef struct {
30: NEP_EXT_OP extop;
31: Vec w;
32: } NEP_SLP_MATSHELL;
34: static PetscErrorCode NEPSetUp_SLP(NEP nep)
35: {
36: NEP_SLP *ctx = (NEP_SLP*)nep->data;
37: PetscBool flg;
38: ST st;
40: PetscFunctionBegin;
41: if (nep->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
42: nep->ncv = nep->nev;
43: if (nep->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
44: nep->mpd = nep->nev;
45: 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: if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
47: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
48: PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
49: NEPCheckUnsupported(nep,NEP_FEATURE_REGION);
51: if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
52: PetscCall(EPSGetST(ctx->eps,&st));
53: PetscCall(PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,""));
54: PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"SLP does not support spectral transformation");
55: PetscCall(EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE));
56: PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_LARGEST_MAGNITUDE));
57: PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(nep->tol)/10.0,nep->max_it));
58: if (nep->tol==(PetscReal)PETSC_DETERMINE) nep->tol = SLEPC_DEFAULT_TOL;
59: if (ctx->deftol==(PetscReal)PETSC_DETERMINE) ctx->deftol = nep->tol;
61: if (nep->twosided) {
62: nep->ops->solve = NEPSolve_SLP_Twosided;
63: nep->ops->computevectors = NULL;
64: if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
65: PetscCall(EPSGetST(ctx->epsts,&st));
66: PetscCall(PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,""));
67: PetscCheck(!flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"SLP does not support spectral transformation");
68: PetscCall(EPSSetDimensions(ctx->epsts,1,PETSC_DECIDE,PETSC_DECIDE));
69: PetscCall(EPSSetWhichEigenpairs(ctx->epsts,EPS_LARGEST_MAGNITUDE));
70: PetscCall(EPSSetTolerances(ctx->epsts,SlepcDefaultTol(nep->tol)/10.0,nep->max_it));
71: } else {
72: nep->ops->solve = NEPSolve_SLP;
73: nep->ops->computevectors = NEPComputeVectors_Schur;
74: }
75: PetscCall(NEPAllocateSolution(nep,0));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode MatMult_SLP(Mat M,Vec x,Vec y)
80: {
81: NEP_SLP_MATSHELL *ctx;
83: PetscFunctionBegin;
84: PetscCall(MatShellGetContext(M,&ctx));
85: PetscCall(MatMult(ctx->extop->MJ,x,ctx->w));
86: PetscCall(NEPDeflationFunctionSolve(ctx->extop,ctx->w,y));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode MatDestroy_SLP(Mat M)
91: {
92: NEP_SLP_MATSHELL *ctx;
94: PetscFunctionBegin;
95: PetscCall(MatShellGetContext(M,&ctx));
96: PetscCall(VecDestroy(&ctx->w));
97: PetscCall(PetscFree(ctx));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
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;
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
114: static PetscErrorCode NEPSLPSetUpLinearEP(NEP nep,NEP_EXT_OP extop,PetscScalar lambda,Vec u,PetscBool ini)
115: {
116: NEP_SLP *slpctx = (NEP_SLP*)nep->data;
117: Mat Mshell;
118: PetscInt nloc,mloc;
119: NEP_SLP_MATSHELL *shellctx;
121: PetscFunctionBegin;
122: if (ini) {
123: /* Create mat shell */
124: PetscCall(PetscNew(&shellctx));
125: shellctx->extop = extop;
126: PetscCall(NEPDeflationCreateVec(extop,&shellctx->w));
127: PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
128: nloc += extop->szd; mloc += extop->szd;
129: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell));
130: PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(PetscErrorCodeFn*)MatMult_SLP));
131: PetscCall(MatShellSetOperation(Mshell,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_SLP));
132: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
133: PetscCall(MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(PetscErrorCodeFn*)MatCreateVecs_SLP));
134: #endif
135: PetscCall(EPSSetOperators(slpctx->eps,Mshell,NULL));
136: PetscCall(MatDestroy(&Mshell));
137: }
138: PetscCall(NEPDeflationSolveSetUp(extop,lambda));
139: PetscCall(NEPDeflationComputeJacobian(extop,lambda,NULL));
140: PetscCall(EPSSetInitialSpace(slpctx->eps,1,&u));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: PetscErrorCode NEPSolve_SLP(NEP nep)
145: {
146: NEP_SLP *ctx = (NEP_SLP*)nep->data;
147: Mat F,H,A;
148: Vec uu,u,r;
149: PetscScalar sigma,lambda,mu,im;
150: PetscReal resnorm;
151: PetscInt nconv;
152: PetscBool skip=PETSC_FALSE,lock=PETSC_FALSE;
153: NEP_EXT_OP extop=NULL; /* Extended operator for deflation */
155: PetscFunctionBegin;
156: /* get initial approximation of eigenvalue and eigenvector */
157: PetscCall(NEPGetDefaultShift(nep,&sigma));
158: if (!nep->nini) PetscCall(BVSetRandomColumn(nep->V,0));
159: lambda = sigma;
160: if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
161: PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_TRUE,nep->nev,&extop));
162: PetscCall(NEPDeflationCreateVec(extop,&u));
163: PetscCall(VecDuplicate(u,&r));
164: PetscCall(BVGetColumn(nep->V,0,&uu));
165: PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE));
166: PetscCall(BVRestoreColumn(nep->V,0,&uu));
168: /* Restart loop */
169: while (nep->reason == NEP_CONVERGED_ITERATING) {
170: nep->its++;
172: /* form residual, r = T(lambda)*u (used in convergence test only) */
173: PetscCall(NEPDeflationComputeFunction(extop,lambda,&F));
174: PetscCall(MatMult(F,u,r));
176: /* convergence test */
177: PetscCall(VecNorm(r,NORM_2,&resnorm));
178: PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
179: nep->eigr[nep->nconv] = lambda;
180: if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
181: if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
182: PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
183: PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
184: PetscCall(MatMult(F,u,r));
185: PetscCall(VecNorm(r,NORM_2,&resnorm));
186: PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
187: if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
188: } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
189: }
191: if (lock) {
192: PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
193: nep->nconv = nep->nconv + 1;
194: skip = PETSC_TRUE;
195: lock = PETSC_FALSE;
196: PetscCall(NEPDeflationLocking(extop,u,lambda));
197: }
198: PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
199: 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));
201: if (nep->reason == NEP_CONVERGED_ITERATING) {
202: if (!skip) {
203: /* evaluate T(lambda) and T'(lambda) */
204: PetscCall(NEPSLPSetUpLinearEP(nep,extop,lambda,u,nep->its==1?PETSC_TRUE:PETSC_FALSE));
205: /* compute new eigenvalue correction mu and eigenvector approximation u */
206: PetscCall(EPSSolve(ctx->eps));
207: PetscCall(EPSGetConverged(ctx->eps,&nconv));
208: if (!nconv) {
209: PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
210: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
211: break;
212: }
213: PetscCall(EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL));
214: mu = 1.0/mu;
215: PetscCheck(PetscAbsScalar(im)<PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
216: } else {
217: nep->its--; /* do not count this as a full iteration */
218: /* use second eigenpair computed in previous iteration */
219: PetscCall(EPSGetConverged(ctx->eps,&nconv));
220: if (nconv>=2) {
221: PetscCall(EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL));
222: mu = 1.0/mu;
223: } else {
224: PetscCall(NEPDeflationSetRandomVec(extop,u));
225: mu = lambda-sigma;
226: }
227: skip = PETSC_FALSE;
228: }
229: /* correct eigenvalue */
230: lambda = lambda - mu;
231: }
232: }
233: PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
234: PetscCall(DSSetType(nep->ds,DSNHEP));
235: PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
236: PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
237: PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
238: PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
239: PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
240: PetscCall(MatDestroy(&H));
241: PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
242: PetscCall(NEPDeflationReset(extop));
243: PetscCall(VecDestroy(&u));
244: PetscCall(VecDestroy(&r));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode NEPSetFromOptions_SLP(NEP nep,PetscOptionItems PetscOptionsObject)
249: {
250: NEP_SLP *ctx = (NEP_SLP*)nep->data;
251: PetscBool flg;
252: PetscReal r;
254: PetscFunctionBegin;
255: PetscOptionsHeadBegin(PetscOptionsObject,"NEP SLP Options");
257: r = 0.0;
258: PetscCall(PetscOptionsReal("-nep_slp_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPSLPSetDeflationThreshold",ctx->deftol,&r,&flg));
259: if (flg) PetscCall(NEPSLPSetDeflationThreshold(nep,r));
261: PetscOptionsHeadEnd();
263: if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
264: PetscCall(EPSSetFromOptions(ctx->eps));
265: if (nep->twosided) {
266: if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
267: PetscCall(EPSSetFromOptions(ctx->epsts));
268: }
269: if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
270: PetscCall(KSPSetFromOptions(ctx->ksp));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode NEPSLPSetDeflationThreshold_SLP(NEP nep,PetscReal deftol)
275: {
276: NEP_SLP *ctx = (NEP_SLP*)nep->data;
278: PetscFunctionBegin;
279: if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) {
280: ctx->deftol = PETSC_DETERMINE;
281: nep->state = NEP_STATE_INITIAL;
282: } else {
283: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
284: ctx->deftol = deftol;
285: }
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: NEPSLPSetDeflationThreshold - Sets the threshold value used to switch between
291: deflated and non-deflated iteration.
293: Logically Collective
295: Input Parameters:
296: + nep - the nonlinear eigensolver context
297: - deftol - the threshold value
299: Options Database Key:
300: . -nep_slp_deflation_threshold \<deftol\> - set the threshold
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.
310: Level: advanced
312: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPGetDeflationThreshold()`
313: @*/
314: PetscErrorCode NEPSLPSetDeflationThreshold(NEP nep,PetscReal deftol)
315: {
316: PetscFunctionBegin;
319: PetscTryMethod(nep,"NEPSLPSetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode NEPSLPGetDeflationThreshold_SLP(NEP nep,PetscReal *deftol)
324: {
325: NEP_SLP *ctx = (NEP_SLP*)nep->data;
327: PetscFunctionBegin;
328: *deftol = ctx->deftol;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*@
333: NEPSLPGetDeflationThreshold - Returns the threshold value that controls deflation.
335: Not Collective
337: Input Parameter:
338: . nep - the nonlinear eigensolver context
340: Output Parameter:
341: . deftol - the threshold
343: Level: advanced
345: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPSetDeflationThreshold()`
346: @*/
347: PetscErrorCode NEPSLPGetDeflationThreshold(NEP nep,PetscReal *deftol)
348: {
349: PetscFunctionBegin;
351: PetscAssertPointer(deftol,2);
352: PetscUseMethod(nep,"NEPSLPGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode NEPSLPSetEPS_SLP(NEP nep,EPS eps)
357: {
358: NEP_SLP *ctx = (NEP_SLP*)nep->data;
360: PetscFunctionBegin;
361: PetscCall(PetscObjectReference((PetscObject)eps));
362: PetscCall(EPSDestroy(&ctx->eps));
363: ctx->eps = eps;
364: nep->state = NEP_STATE_INITIAL;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@
369: NEPSLPSetEPS - Associate a linear eigensolver object (`EPS`) to the
370: nonlinear eigenvalue solver.
372: Collective
374: Input Parameters:
375: + nep - the nonlinear eigensolver context
376: - eps - the linear eigensolver context
378: Level: advanced
380: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPGetEPS()`
381: @*/
382: PetscErrorCode NEPSLPSetEPS(NEP nep,EPS eps)
383: {
384: PetscFunctionBegin;
387: PetscCheckSameComm(nep,1,eps,2);
388: PetscTryMethod(nep,"NEPSLPSetEPS_C",(NEP,EPS),(nep,eps));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode NEPSLPGetEPS_SLP(NEP nep,EPS *eps)
393: {
394: NEP_SLP *ctx = (NEP_SLP*)nep->data;
396: PetscFunctionBegin;
397: if (!ctx->eps) {
398: PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps));
399: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1));
400: PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix));
401: PetscCall(EPSAppendOptionsPrefix(ctx->eps,"nep_slp_"));
402: PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options));
403: }
404: *eps = ctx->eps;
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: NEPSLPGetEPS - Retrieve the linear eigensolver object (`EPS`) associated
410: to the nonlinear eigenvalue solver.
412: Collective
414: Input Parameter:
415: . nep - the nonlinear eigensolver context
417: Output Parameter:
418: . eps - the linear eigensolver context
420: Level: advanced
422: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPSetEPS()`
423: @*/
424: PetscErrorCode NEPSLPGetEPS(NEP nep,EPS *eps)
425: {
426: PetscFunctionBegin;
428: PetscAssertPointer(eps,2);
429: PetscUseMethod(nep,"NEPSLPGetEPS_C",(NEP,EPS*),(nep,eps));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static PetscErrorCode NEPSLPSetEPSLeft_SLP(NEP nep,EPS eps)
434: {
435: NEP_SLP *ctx = (NEP_SLP*)nep->data;
437: PetscFunctionBegin;
438: PetscCall(PetscObjectReference((PetscObject)eps));
439: PetscCall(EPSDestroy(&ctx->epsts));
440: ctx->epsts = eps;
441: nep->state = NEP_STATE_INITIAL;
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
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.
450: Collective
452: Input Parameters:
453: + nep - the nonlinear eigensolver context
454: - eps - the linear eigensolver context
456: Level: advanced
458: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPGetEPSLeft()`, `NEPSetTwoSided()`
459: @*/
460: PetscErrorCode NEPSLPSetEPSLeft(NEP nep,EPS eps)
461: {
462: PetscFunctionBegin;
465: PetscCheckSameComm(nep,1,eps,2);
466: PetscTryMethod(nep,"NEPSLPSetEPSLeft_C",(NEP,EPS),(nep,eps));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: static PetscErrorCode NEPSLPGetEPSLeft_SLP(NEP nep,EPS *eps)
471: {
472: NEP_SLP *ctx = (NEP_SLP*)nep->data;
474: PetscFunctionBegin;
475: if (!ctx->epsts) {
476: PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->epsts));
477: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->epsts,(PetscObject)nep,1));
478: PetscCall(EPSSetOptionsPrefix(ctx->epsts,((PetscObject)nep)->prefix));
479: PetscCall(EPSAppendOptionsPrefix(ctx->epsts,"nep_slp_left_"));
480: PetscCall(PetscObjectSetOptions((PetscObject)ctx->epsts,((PetscObject)nep)->options));
481: }
482: *eps = ctx->epsts;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
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.
491: Collective
493: Input Parameter:
494: . nep - the nonlinear eigensolver context
496: Output Parameter:
497: . eps - the linear eigensolver context
499: Level: advanced
501: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPSetEPSLeft()`, `NEPSetTwoSided()`
502: @*/
503: PetscErrorCode NEPSLPGetEPSLeft(NEP nep,EPS *eps)
504: {
505: PetscFunctionBegin;
507: PetscAssertPointer(eps,2);
508: PetscUseMethod(nep,"NEPSLPGetEPSLeft_C",(NEP,EPS*),(nep,eps));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode NEPSLPSetKSP_SLP(NEP nep,KSP ksp)
513: {
514: NEP_SLP *ctx = (NEP_SLP*)nep->data;
516: PetscFunctionBegin;
517: PetscCall(PetscObjectReference((PetscObject)ksp));
518: PetscCall(KSPDestroy(&ctx->ksp));
519: ctx->ksp = ksp;
520: nep->state = NEP_STATE_INITIAL;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: NEPSLPSetKSP - Associate a linear solver object (`KSP`) to the nonlinear
526: eigenvalue solver.
528: Collective
530: Input Parameters:
531: + nep - the nonlinear eigensolver context
532: - ksp - the linear solver object
534: Note:
535: This `KSP` object is used only for deflation, when computing more that
536: one eigenpair.
538: Level: advanced
540: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPGetKSP()`
541: @*/
542: PetscErrorCode NEPSLPSetKSP(NEP nep,KSP ksp)
543: {
544: PetscFunctionBegin;
547: PetscCheckSameComm(nep,1,ksp,2);
548: PetscTryMethod(nep,"NEPSLPSetKSP_C",(NEP,KSP),(nep,ksp));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
553: {
554: NEP_SLP *ctx = (NEP_SLP*)nep->data;
556: PetscFunctionBegin;
557: if (!ctx->ksp) {
558: PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
559: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
560: PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
561: PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_"));
562: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
563: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
564: PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
565: }
566: *ksp = ctx->ksp;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@
571: NEPSLPGetKSP - Retrieve the linear solver object (`KSP`) associated with
572: the nonlinear eigenvalue solver.
574: Collective
576: Input Parameter:
577: . nep - the nonlinear eigensolver context
579: Output Parameter:
580: . ksp - the linear solver object
582: Note:
583: This `KSP` object is used only for deflation, when computing more that
584: one eigenpair.
586: Level: advanced
588: .seealso: [](ch:nep), `NEPSLP`, `NEPSLPSetKSP()`
589: @*/
590: PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
591: {
592: PetscFunctionBegin;
594: PetscAssertPointer(ksp,2);
595: PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode NEPView_SLP(NEP nep,PetscViewer viewer)
600: {
601: NEP_SLP *ctx = (NEP_SLP*)nep->data;
602: PetscBool isascii;
604: PetscFunctionBegin;
605: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
606: if (isascii) {
607: if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer," deflation threshold: %g\n",(double)ctx->deftol));
608: if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
609: PetscCall(PetscViewerASCIIPushTab(viewer));
610: PetscCall(EPSView(ctx->eps,viewer));
611: if (nep->twosided) {
612: if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
613: PetscCall(EPSView(ctx->epsts,viewer));
614: }
615: if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
616: PetscCall(KSPView(ctx->ksp,viewer));
617: PetscCall(PetscViewerASCIIPopTab(viewer));
618: }
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: static PetscErrorCode NEPReset_SLP(NEP nep)
623: {
624: NEP_SLP *ctx = (NEP_SLP*)nep->data;
626: PetscFunctionBegin;
627: PetscCall(EPSReset(ctx->eps));
628: if (nep->twosided) PetscCall(EPSReset(ctx->epsts));
629: PetscCall(KSPReset(ctx->ksp));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: static PetscErrorCode NEPDestroy_SLP(NEP nep)
634: {
635: NEP_SLP *ctx = (NEP_SLP*)nep->data;
637: PetscFunctionBegin;
638: PetscCall(KSPDestroy(&ctx->ksp));
639: PetscCall(EPSDestroy(&ctx->eps));
640: PetscCall(EPSDestroy(&ctx->epsts));
641: PetscCall(PetscFree(nep->data));
642: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NULL));
643: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NULL));
644: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL));
645: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL));
646: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NULL));
647: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NULL));
648: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL));
649: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: /*MC
654: NEPSLP - NEPSLP = "slp" - Successive Linear Problems method.
656: Notes:
657: This solver is based on the classical Successive Linear Problems
658: (SLP) method proposed by {cite:t}`Ruh73`. At each step, this
659: method has to solve a linear eigenvalue problem. Call
660: `NEPSLPGetEPS()` to configure the `EPS` object used for this.
662: `NEPSLP` supports computing left eigenvectors when `NEPSetTwoSided()`
663: has been set. In that case, a different `EPS` is used for the left
664: recurrence, see `NEPSLPGetEPSLeft()`.
666: The solver incorporates deflation, so that several eigenpairs con be
667: computed. Details of the implementation in SLEPc can be found in
668: {cite:p}`Cam21`.
670: Level: beginner
672: .seealso: [](ch:nep), `NEP`, `NEPType`, `NEPSetType()`, `NEPSLPGetEPS()`, `NEPSetTwoSided()`, `NEPSLPGetEPSLeft()`
673: M*/
674: SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
675: {
676: NEP_SLP *ctx;
678: PetscFunctionBegin;
679: PetscCall(PetscNew(&ctx));
680: nep->data = (void*)ctx;
682: nep->useds = PETSC_TRUE;
683: ctx->deftol = PETSC_DETERMINE;
685: nep->ops->solve = NEPSolve_SLP;
686: nep->ops->setup = NEPSetUp_SLP;
687: nep->ops->setfromoptions = NEPSetFromOptions_SLP;
688: nep->ops->reset = NEPReset_SLP;
689: nep->ops->destroy = NEPDestroy_SLP;
690: nep->ops->view = NEPView_SLP;
691: nep->ops->computevectors = NEPComputeVectors_Schur;
693: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NEPSLPSetDeflationThreshold_SLP));
694: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NEPSLPGetDeflationThreshold_SLP));
695: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP));
696: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP));
697: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NEPSLPSetEPSLeft_SLP));
698: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NEPSLPGetEPSLeft_SLP));
699: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP));
700: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }