Actual source code: slp.c

slepc-3.20.1 2023-11-27
Report Typos and Errors
  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_DEFAULT) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
 42:   nep->ncv = nep->nev;
 43:   if (nep->mpd!=PETSC_DEFAULT) 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_DEFAULT) 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_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
 59:   if (ctx->deftol==(PetscReal)PETSC_DEFAULT) 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)
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,(void(*)(void))MatMult_SLP));
131:     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:     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) {
280:     ctx->deftol = PETSC_DEFAULT;
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    - nonlinear eigenvalue solver
297: -  deftol - the threshold value

299:    Options Database Keys:
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: 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 - nonlinear eigenvalue solver

340:    Output Parameter:
341: .  deftol - the threshold

343:    Level: advanced

345: .seealso: 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 - nonlinear eigenvalue solver
376: -  eps - the eigensolver object

378:    Level: advanced

380: .seealso: 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 - nonlinear eigenvalue solver

417:    Output Parameter:
418: .  eps - the eigensolver object

420:    Level: advanced

422: .seealso: 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 - nonlinear eigenvalue solver
454: -  eps - the eigensolver object

456:    Level: advanced

458: .seealso: 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 - nonlinear eigenvalue solver

496:    Output Parameter:
497: .  eps - the eigensolver object

499:    Level: advanced

501: .seealso: 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 - eigenvalue solver
532: -  ksp - the linear solver object

534:    Level: advanced

536: .seealso: NEPSLPGetKSP()
537: @*/
538: PetscErrorCode NEPSLPSetKSP(NEP nep,KSP ksp)
539: {
540:   PetscFunctionBegin;
543:   PetscCheckSameComm(nep,1,ksp,2);
544:   PetscTryMethod(nep,"NEPSLPSetKSP_C",(NEP,KSP),(nep,ksp));
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
549: {
550:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

552:   PetscFunctionBegin;
553:   if (!ctx->ksp) {
554:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
555:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
556:     PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
557:     PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_"));
558:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
559:     PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
560:     PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
561:   }
562:   *ksp = ctx->ksp;
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: /*@
567:    NEPSLPGetKSP - Retrieve the linear solver object (KSP) associated with
568:    the nonlinear eigenvalue solver.

570:    Collective

572:    Input Parameter:
573: .  nep - nonlinear eigenvalue solver

575:    Output Parameter:
576: .  ksp - the linear solver object

578:    Level: advanced

580: .seealso: NEPSLPSetKSP()
581: @*/
582: PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
583: {
584:   PetscFunctionBegin;
586:   PetscAssertPointer(ksp,2);
587:   PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: static PetscErrorCode NEPView_SLP(NEP nep,PetscViewer viewer)
592: {
593:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
594:   PetscBool      isascii;

596:   PetscFunctionBegin;
597:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
598:   if (isascii) {
599:     if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer,"  deflation threshold: %g\n",(double)ctx->deftol));
600:     if (!ctx->eps) PetscCall(NEPSLPGetEPS(nep,&ctx->eps));
601:     PetscCall(PetscViewerASCIIPushTab(viewer));
602:     PetscCall(EPSView(ctx->eps,viewer));
603:     if (nep->twosided) {
604:       if (!ctx->epsts) PetscCall(NEPSLPGetEPSLeft(nep,&ctx->epsts));
605:       PetscCall(EPSView(ctx->epsts,viewer));
606:     }
607:     if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
608:     PetscCall(KSPView(ctx->ksp,viewer));
609:     PetscCall(PetscViewerASCIIPopTab(viewer));
610:   }
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: static PetscErrorCode NEPReset_SLP(NEP nep)
615: {
616:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

618:   PetscFunctionBegin;
619:   PetscCall(EPSReset(ctx->eps));
620:   if (nep->twosided) PetscCall(EPSReset(ctx->epsts));
621:   PetscCall(KSPReset(ctx->ksp));
622:   PetscFunctionReturn(PETSC_SUCCESS);
623: }

625: static PetscErrorCode NEPDestroy_SLP(NEP nep)
626: {
627:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

629:   PetscFunctionBegin;
630:   PetscCall(KSPDestroy(&ctx->ksp));
631:   PetscCall(EPSDestroy(&ctx->eps));
632:   PetscCall(EPSDestroy(&ctx->epsts));
633:   PetscCall(PetscFree(nep->data));
634:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NULL));
635:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NULL));
636:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL));
637:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL));
638:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NULL));
639:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NULL));
640:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL));
641:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
646: {
647:   NEP_SLP        *ctx;

649:   PetscFunctionBegin;
650:   PetscCall(PetscNew(&ctx));
651:   nep->data = (void*)ctx;

653:   nep->useds  = PETSC_TRUE;
654:   ctx->deftol = PETSC_DEFAULT;

656:   nep->ops->solve          = NEPSolve_SLP;
657:   nep->ops->setup          = NEPSetUp_SLP;
658:   nep->ops->setfromoptions = NEPSetFromOptions_SLP;
659:   nep->ops->reset          = NEPReset_SLP;
660:   nep->ops->destroy        = NEPDestroy_SLP;
661:   nep->ops->view           = NEPView_SLP;
662:   nep->ops->computevectors = NEPComputeVectors_Schur;

664:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetDeflationThreshold_C",NEPSLPSetDeflationThreshold_SLP));
665:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetDeflationThreshold_C",NEPSLPGetDeflationThreshold_SLP));
666:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP));
667:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP));
668:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPSLeft_C",NEPSLPSetEPSLeft_SLP));
669:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPSLeft_C",NEPSLPGetEPSLeft_SLP));
670:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP));
671:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP));
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }