Actual source code: nleigs-fullb.c

slepc-3.22.2 2024-12-02
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:    Full basis for the linearization of the rational approximation of non-linear eigenproblems
 12: */

 14: #include <slepc/private/nepimpl.h>
 15: #include "nleigs.h"

 17: static PetscErrorCode MatMult_FullBasis_Sinvert(Mat M,Vec x,Vec y)
 18: {
 19:   NEP_NLEIGS        *ctx;
 20:   NEP               nep;
 21:   const PetscScalar *px;
 22:   PetscScalar       *beta,*s,*xi,*t,*py,sigma;
 23:   PetscInt          nmat,d,i,k,m;
 24:   Vec               xx,xxx,yy,yyy,w,ww,www;

 26:   PetscFunctionBegin;
 27:   PetscCall(MatShellGetContext(M,&nep));
 28:   ctx = (NEP_NLEIGS*)nep->data;
 29:   beta = ctx->beta; s = ctx->s; xi = ctx->xi;
 30:   sigma = ctx->shifts[0];
 31:   nmat = ctx->nmat;
 32:   d = nmat-1;
 33:   m = nep->nloc;
 34:   PetscCall(PetscMalloc1(ctx->nmat,&t));
 35:   xx = ctx->w[0]; xxx = ctx->w[1]; yy = ctx->w[2]; yyy=ctx->w[3];
 36:   w = nep->work[0]; ww = nep->work[1]; www = nep->work[2];
 37:   PetscCall(VecGetArrayRead(x,&px));
 38:   PetscCall(VecGetArray(y,&py));
 39:   PetscCall(VecPlaceArray(xx,px+(d-1)*m));
 40:   PetscCall(VecPlaceArray(xxx,px+(d-2)*m));
 41:   PetscCall(VecPlaceArray(yy,py+(d-2)*m));
 42:   PetscCall(VecCopy(xxx,yy));
 43:   PetscCall(VecAXPY(yy,beta[d-1]/xi[d-2],xx));
 44:   PetscCall(VecScale(yy,1.0/(s[d-2]-sigma)));
 45:   PetscCall(VecResetArray(xx));
 46:   PetscCall(VecResetArray(xxx));
 47:   PetscCall(VecResetArray(yy));
 48:   for (i=d-3;i>=0;i--) {
 49:     PetscCall(VecPlaceArray(xx,px+(i+1)*m));
 50:     PetscCall(VecPlaceArray(xxx,px+i*m));
 51:     PetscCall(VecPlaceArray(yy,py+i*m));
 52:     PetscCall(VecPlaceArray(yyy,py+(i+1)*m));
 53:     PetscCall(VecCopy(xxx,yy));
 54:     PetscCall(VecAXPY(yy,beta[i+1]/xi[i],xx));
 55:     PetscCall(VecAXPY(yy,-beta[i+1]*(1.0-sigma/xi[i]),yyy));
 56:     PetscCall(VecScale(yy,1.0/(s[i]-sigma)));
 57:     PetscCall(VecResetArray(xx));
 58:     PetscCall(VecResetArray(xxx));
 59:     PetscCall(VecResetArray(yy));
 60:     PetscCall(VecResetArray(yyy));
 61:   }
 62:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
 63:     PetscCall(VecZeroEntries(w));
 64:     for (k=0;k<nep->nt;k++) {
 65:       PetscCall(VecZeroEntries(ww));
 66:       PetscCall(VecPlaceArray(xx,px+(d-1)*m));
 67:       PetscCall(VecAXPY(ww,-ctx->coeffD[k+nep->nt*d]/beta[d],xx));
 68:       PetscCall(VecResetArray(xx));
 69:       for (i=0;i<d-1;i++) {
 70:         PetscCall(VecPlaceArray(yy,py+i*m));
 71:         PetscCall(VecAXPY(ww,-ctx->coeffD[nep->nt*i+k],yy));
 72:         PetscCall(VecResetArray(yy));
 73:       }
 74:       PetscCall(MatMult(nep->A[k],ww,www));
 75:       PetscCall(VecAXPY(w,1.0,www));
 76:     }
 77:   } else {
 78:     PetscCall(VecPlaceArray(xx,px+(d-1)*m));
 79:     PetscCall(MatMult(ctx->D[d],xx,w));
 80:     PetscCall(VecScale(w,-1.0/beta[d]));
 81:     PetscCall(VecResetArray(xx));
 82:     for (i=0;i<d-1;i++) {
 83:       PetscCall(VecPlaceArray(yy,py+i*m));
 84:       PetscCall(MatMult(ctx->D[i],yy,ww));
 85:       PetscCall(VecResetArray(yy));
 86:       PetscCall(VecAXPY(w,-1.0,ww));
 87:     }
 88:   }
 89:   PetscCall(VecPlaceArray(yy,py+(d-1)*m));
 90:   PetscCall(KSPSolve(ctx->ksp[0],w,yy));
 91:   PetscCall(NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t));
 92:   for (i=0;i<d-1;i++) {
 93:     PetscCall(VecPlaceArray(yyy,py+i*m));
 94:     PetscCall(VecAXPY(yyy,t[i],yy));
 95:     PetscCall(VecResetArray(yyy));
 96:   }
 97:   PetscCall(VecScale(yy,t[d-1]));
 98:   PetscCall(VecResetArray(yy));
 99:   PetscCall(VecRestoreArrayRead(x,&px));
100:   PetscCall(VecRestoreArray(y,&py));
101:   PetscCall(PetscFree(t));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: static PetscErrorCode MatMultTranspose_FullBasis_Sinvert(Mat M,Vec x,Vec y)
106: {
107:   NEP_NLEIGS        *ctx;
108:   NEP               nep;
109:   const PetscScalar *px;
110:   PetscScalar       *beta,*s,*xi,*t,*py,sigma;
111:   PetscInt          nmat,d,i,k,m;
112:   Vec               xx,yy,yyy,w,z0;

114:   PetscFunctionBegin;
115:   PetscCall(MatShellGetContext(M,&nep));
116:   ctx = (NEP_NLEIGS*)nep->data;
117:   beta = ctx->beta; s = ctx->s; xi = ctx->xi;
118:   sigma = ctx->shifts[0];
119:   nmat = ctx->nmat;
120:   d = nmat-1;
121:   m = nep->nloc;
122:   PetscCall(PetscMalloc1(ctx->nmat,&t));
123:   xx = ctx->w[0]; yy = ctx->w[1]; yyy=ctx->w[2];
124:   w = nep->work[0]; z0 = nep->work[1];
125:   PetscCall(VecGetArrayRead(x,&px));
126:   PetscCall(VecGetArray(y,&py));
127:   PetscCall(NEPNLEIGSEvalNRTFunct(nep,d,sigma,t));
128:   PetscCall(VecPlaceArray(xx,px+(d-1)*m));
129:   PetscCall(VecCopy(xx,w));
130:   PetscCall(VecScale(w,t[d-1]));
131:   PetscCall(VecResetArray(xx));
132:   for (i=0;i<d-1;i++) {
133:     PetscCall(VecPlaceArray(xx,px+i*m));
134:     PetscCall(VecAXPY(w,t[i],xx));
135:     PetscCall(VecResetArray(xx));
136:   }
137:   PetscCall(KSPSolveTranspose(ctx->ksp[0],w,z0));

139:   PetscCall(VecPlaceArray(yy,py));
140:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
141:     PetscCall(VecZeroEntries(yy));
142:     for (k=0;k<nep->nt;k++) {
143:       PetscCall(MatMult(nep->A[k],z0,w));
144:       PetscCall(VecAXPY(yy,ctx->coeffD[k],w));
145:     }
146:   } else PetscCall(MatMultTranspose(ctx->D[0],z0,yy));
147:   PetscCall(VecPlaceArray(xx,px));
148:   PetscCall(VecAXPY(yy,-1.0,xx));
149:   PetscCall(VecResetArray(xx));
150:   PetscCall(VecScale(yy,-1.0/(s[0]-sigma)));
151:   PetscCall(VecResetArray(yy));
152:   for (i=2;i<d;i++) {
153:     PetscCall(VecPlaceArray(yy,py+(i-1)*m));
154:     if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
155:       PetscCall(VecZeroEntries(yy));
156:       for (k=0;k<nep->nt;k++) {
157:         PetscCall(MatMult(nep->A[k],z0,w));
158:         PetscCall(VecAXPY(yy,ctx->coeffD[k+(i-1)*nep->nt],w));
159:       }
160:     } else PetscCall(MatMultTranspose(ctx->D[i-1],z0,yy));
161:     PetscCall(VecPlaceArray(yyy,py+(i-2)*m));
162:     PetscCall(VecAXPY(yy,beta[i-1]*(1.0-sigma/xi[i-2]),yyy));
163:     PetscCall(VecResetArray(yyy));
164:     PetscCall(VecPlaceArray(xx,px+(i-1)*m));
165:     PetscCall(VecAXPY(yy,-1.0,xx));
166:     PetscCall(VecResetArray(xx));
167:     PetscCall(VecScale(yy,-1.0/(s[i-1]-sigma)));
168:     PetscCall(VecResetArray(yy));
169:   }
170:   PetscCall(VecPlaceArray(yy,py+(d-1)*m));
171:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
172:     PetscCall(VecZeroEntries(yy));
173:     for (k=0;k<nep->nt;k++) {
174:       PetscCall(MatMult(nep->A[k],z0,w));
175:       PetscCall(VecAXPY(yy,ctx->coeffD[k+d*nep->nt],w));
176:     }
177:   } else PetscCall(MatMultTranspose(ctx->D[d],z0,yy));
178:   PetscCall(VecScale(yy,-1.0/beta[d]));
179:   PetscCall(VecPlaceArray(yyy,py+(d-2)*m));
180:   PetscCall(VecAXPY(yy,beta[d-1]/xi[d-2],yyy));
181:   PetscCall(VecResetArray(yyy));
182:   PetscCall(VecResetArray(yy));

184:   for (i=d-2;i>0;i--) {
185:     PetscCall(VecPlaceArray(yyy,py+(i-1)*m));
186:     PetscCall(VecPlaceArray(yy,py+i*m));
187:     PetscCall(VecAXPY(yy,beta[i]/xi[i-1],yyy));
188:     PetscCall(VecResetArray(yyy));
189:     PetscCall(VecResetArray(yy));
190:   }

192:   PetscCall(VecRestoreArrayRead(x,&px));
193:   PetscCall(VecRestoreArray(y,&py));
194:   PetscCall(PetscFree(t));
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: static PetscErrorCode BackTransform_FullBasis(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
199: {
200:   NEP            nep;

202:   PetscFunctionBegin;
203:   PetscCall(STShellGetContext(st,&nep));
204:   PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,n,eigr,eigi));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode Apply_FullBasis(ST st,Vec x,Vec y)
209: {
210:   NEP            nep;
211:   NEP_NLEIGS     *ctx;

213:   PetscFunctionBegin;
214:   PetscCall(STShellGetContext(st,&nep));
215:   ctx = (NEP_NLEIGS*)nep->data;
216:   PetscCall(MatMult(ctx->A,x,y));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: static PetscErrorCode ApplyTranspose_FullBasis(ST st,Vec x,Vec y)
221: {
222:   NEP            nep;
223:   NEP_NLEIGS     *ctx;

225:   PetscFunctionBegin;
226:   PetscCall(STShellGetContext(st,&nep));
227:   ctx = (NEP_NLEIGS*)nep->data;
228:   PetscCall(MatMultTranspose(ctx->A,x,y));
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: PetscErrorCode NEPSetUp_NLEIGS_FullBasis(NEP nep)
233: {
234:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
235:   ST             st;
236:   Mat            Q;
237:   PetscInt       i=0,deg=ctx->nmat-1;
238:   PetscBool      trackall,istrivial,ks;
239:   PetscScalar    *epsarray,*neparray;
240:   Vec            veps,w=NULL;
241:   EPSWhich       which;

243:   PetscFunctionBegin;
244:   PetscCheck(ctx->nshifts==0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The full-basis option is not supported with rational Krylov");
245:   if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
246:   PetscCall(EPSGetST(ctx->eps,&st));
247:   PetscCall(EPSSetTarget(ctx->eps,nep->target));
248:   PetscCall(STSetDefaultShift(st,nep->target));
249:   if (!((PetscObject)ctx->eps)->type_name) PetscCall(EPSSetType(ctx->eps,EPSKRYLOVSCHUR));
250:   else {
251:     PetscCall(PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks));
252:     PetscCheck(ks,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Full-basis option only implemented for Krylov-Schur");
253:   }
254:   PetscCall(STSetType(st,STSHELL));
255:   PetscCall(STShellSetContext(st,nep));
256:   PetscCall(STShellSetBackTransform(st,BackTransform_FullBasis));
257:   PetscCall(KSPGetOperators(ctx->ksp[0],&Q,NULL));
258:   PetscCall(MatCreateVecsEmpty(Q,&ctx->w[0],&ctx->w[1]));
259:   PetscCall(MatCreateVecsEmpty(Q,&ctx->w[2],&ctx->w[3]));
260:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),deg*nep->nloc,deg*nep->nloc,deg*nep->n,deg*nep->n,nep,&ctx->A));
261:   PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_FullBasis_Sinvert));
262:   PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_FullBasis_Sinvert));
263:   PetscCall(STShellSetApply(st,Apply_FullBasis));
264:   PetscCall(STShellSetApplyTranspose(st,ApplyTranspose_FullBasis));
265:   PetscCall(EPSSetOperators(ctx->eps,ctx->A,NULL));
266:   PetscCall(EPSSetProblemType(ctx->eps,EPS_NHEP));
267:   switch (nep->which) {
268:     case NEP_TARGET_MAGNITUDE:   which = EPS_TARGET_MAGNITUDE; break;
269:     case NEP_TARGET_REAL:        which = EPS_TARGET_REAL; break;
270:     case NEP_TARGET_IMAGINARY:   which = EPS_TARGET_IMAGINARY; break;
271:     case NEP_WHICH_USER:         which = EPS_WHICH_USER;
272:       PetscCall(EPSSetEigenvalueComparison(ctx->eps,nep->sc->comparison,nep->sc->comparisonctx));
273:       break;
274:     default: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Should set a target selection in NEPSetWhichEigenpairs()");
275:   }
276:   PetscCall(EPSSetWhichEigenpairs(ctx->eps,which));
277:   PetscCall(RGIsTrivial(nep->rg,&istrivial));
278:   if (!istrivial) PetscCall(EPSSetRG(ctx->eps,nep->rg));
279:   PetscCall(EPSSetDimensions(ctx->eps,nep->nev,nep->ncv,nep->mpd));
280:   PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(nep->tol),nep->max_it));
281:   PetscCall(EPSSetTwoSided(ctx->eps,nep->twosided));
282:   /* Transfer the trackall option from pep to eps */
283:   PetscCall(NEPGetTrackAll(nep,&trackall));
284:   PetscCall(EPSSetTrackAll(ctx->eps,trackall));

286:   /* process initial vector */
287:   if (nep->nini<0) {
288:     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*nep->nloc,deg*nep->n,&veps));
289:     PetscCall(VecGetArray(veps,&epsarray));
290:     for (i=0;i<deg;i++) {
291:       if (i<-nep->nini) {
292:         PetscCall(VecGetArray(nep->IS[i],&neparray));
293:         PetscCall(PetscArraycpy(epsarray+i*nep->nloc,neparray,nep->nloc));
294:         PetscCall(VecRestoreArray(nep->IS[i],&neparray));
295:       } else {
296:         if (!w) PetscCall(VecDuplicate(nep->IS[0],&w));
297:         PetscCall(VecSetRandom(w,NULL));
298:         PetscCall(VecGetArray(w,&neparray));
299:         PetscCall(PetscArraycpy(epsarray+i*nep->nloc,neparray,nep->nloc));
300:         PetscCall(VecRestoreArray(w,&neparray));
301:       }
302:     }
303:     PetscCall(VecRestoreArray(veps,&epsarray));
304:     PetscCall(EPSSetInitialSpace(ctx->eps,1,&veps));
305:     PetscCall(VecDestroy(&veps));
306:     PetscCall(VecDestroy(&w));
307:     PetscCall(SlepcBasisDestroy_Private(&nep->nini,&nep->IS));
308:   }

310:   PetscCall(EPSSetUp(ctx->eps));
311:   PetscCall(EPSGetDimensions(ctx->eps,NULL,&nep->ncv,&nep->mpd));
312:   PetscCall(EPSGetTolerances(ctx->eps,NULL,&nep->max_it));
313:   PetscCall(NEPAllocateSolution(nep,0));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: /*
318:    NEPNLEIGSExtract_None - Extracts the first block of the basis
319:    and normalizes the columns.
320: */
321: static PetscErrorCode NEPNLEIGSExtract_None(NEP nep,EPS eps)
322: {
323:   PetscInt          i,k,m,d;
324:   const PetscScalar *px;
325:   PetscScalar       sigma=nep->target,*b;
326:   Mat               A;
327:   Vec               xxr,xxi=NULL,w,t,xx;
328:   PetscReal         norm;
329:   NEP_NLEIGS        *ctx=(NEP_NLEIGS*)nep->data;

331:   PetscFunctionBegin;
332:   d = ctx->nmat-1;
333:   PetscCall(EPSGetOperators(eps,&A,NULL));
334:   PetscCall(MatCreateVecs(A,&xxr,NULL));
335: #if !defined(PETSC_USE_COMPLEX)
336:   PetscCall(VecDuplicate(xxr,&xxi));
337: #endif
338:   w = nep->work[0];
339:   for (i=0;i<nep->nconv;i++) {
340:     PetscCall(EPSGetEigenvector(eps,i,xxr,xxi));
341:     PetscCall(VecGetArrayRead(xxr,&px));
342:     PetscCall(VecPlaceArray(w,px));
343:     PetscCall(BVInsertVec(nep->V,i,w));
344:     PetscCall(BVNormColumn(nep->V,i,NORM_2,&norm));
345:     PetscCall(BVScaleColumn(nep->V,i,1.0/norm));
346:     PetscCall(VecResetArray(w));
347:     PetscCall(VecRestoreArrayRead(xxr,&px));
348:   }
349:   if (nep->twosided) {
350:     PetscCall(PetscMalloc1(ctx->nmat,&b));
351:     PetscCall(NEPNLEIGSEvalNRTFunct(nep,d,sigma,b));
352:     m = nep->nloc;
353:     xx = ctx->w[0];
354:     w = nep->work[0]; t = nep->work[1];
355:     for (k=0;k<nep->nconv;k++) {
356:       PetscCall(EPSGetLeftEigenvector(eps,k,xxr,xxi));
357:       PetscCall(VecGetArrayRead(xxr,&px));
358:       PetscCall(VecPlaceArray(xx,px+(d-1)*m));
359:       PetscCall(VecCopy(xx,w));
360:       PetscCall(VecScale(w,PetscConj(b[d-1])));
361:       PetscCall(VecResetArray(xx));
362:       for (i=0;i<d-1;i++) {
363:         PetscCall(VecPlaceArray(xx,px+i*m));
364:         PetscCall(VecAXPY(w,PetscConj(b[i]),xx));
365:         PetscCall(VecResetArray(xx));
366:       }
367:       PetscCall(VecConjugate(w));
368:       PetscCall(KSPSolveTranspose(ctx->ksp[0],w,t));
369:       PetscCall(VecConjugate(t));
370:       PetscCall(BVInsertVec(nep->W,k,t));
371:       PetscCall(BVNormColumn(nep->W,k,NORM_2,&norm));
372:       PetscCall(BVScaleColumn(nep->W,k,1.0/norm));
373:       PetscCall(VecRestoreArrayRead(xxr,&px));
374:     }
375:     PetscCall(PetscFree(b));
376:   }
377:   PetscCall(VecDestroy(&xxr));
378: #if !defined(PETSC_USE_COMPLEX)
379:   PetscCall(VecDestroy(&xxi));
380: #endif
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: PetscErrorCode NEPSolve_NLEIGS_FullBasis(NEP nep)
385: {
386:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
387:   PetscInt       i;
388:   PetscScalar    eigi=0.0;

390:   PetscFunctionBegin;
391:   PetscCall(EPSSolve(ctx->eps));
392:   PetscCall(EPSGetConverged(ctx->eps,&nep->nconv));
393:   PetscCall(EPSGetIterationNumber(ctx->eps,&nep->its));
394:   PetscCall(EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&nep->reason));

396:   /* recover eigenvalues */
397:   for (i=0;i<nep->nconv;i++) {
398:     PetscCall(EPSGetEigenpair(ctx->eps,i,&nep->eigr[i],&eigi,NULL,NULL));
399: #if !defined(PETSC_USE_COMPLEX)
400:     PetscCheck(eigi==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex value requires complex arithmetic");
401: #endif
402:   }
403:   PetscCall(NEPNLEIGSExtract_None(nep,ctx->eps));
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

407: PetscErrorCode NEPNLEIGSSetEPS_NLEIGS(NEP nep,EPS eps)
408: {
409:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

411:   PetscFunctionBegin;
412:   PetscCall(PetscObjectReference((PetscObject)eps));
413:   PetscCall(EPSDestroy(&ctx->eps));
414:   ctx->eps = eps;
415:   nep->state = NEP_STATE_INITIAL;
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: /*@
420:    NEPNLEIGSSetEPS - Associate an eigensolver object (EPS) to the NLEIGS solver.

422:    Collective

424:    Input Parameters:
425: +  nep - nonlinear eigenvalue solver
426: -  eps - the eigensolver object

428:    Level: advanced

430: .seealso: NEPNLEIGSGetEPS()
431: @*/
432: PetscErrorCode NEPNLEIGSSetEPS(NEP nep,EPS eps)
433: {
434:   PetscFunctionBegin;
437:   PetscCheckSameComm(nep,1,eps,2);
438:   PetscTryMethod(nep,"NEPNLEIGSSetEPS_C",(NEP,EPS),(nep,eps));
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: static PetscErrorCode EPSMonitor_NLEIGS(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
443: {
444:   NEP            nep = (NEP)ctx;
445:   PetscInt       i,nv = PetscMin(nest,nep->ncv);

447:   PetscFunctionBegin;
448:   for (i=0;i<nv;i++) {
449:     nep->eigr[i]   = eigr[i];
450:     nep->eigi[i]   = eigi[i];
451:     nep->errest[i] = errest[i];
452:   }
453:   PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nv,nep->eigr,nep->eigi));
454:   PetscCall(NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest));
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: PetscErrorCode NEPNLEIGSGetEPS_NLEIGS(NEP nep,EPS *eps)
459: {
460:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

462:   PetscFunctionBegin;
463:   if (!ctx->eps) {
464:     PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps));
465:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1));
466:     PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix));
467:     PetscCall(EPSAppendOptionsPrefix(ctx->eps,"nep_nleigs_"));
468:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options));
469:     PetscCall(EPSMonitorSet(ctx->eps,EPSMonitor_NLEIGS,nep,NULL));
470:   }
471:   *eps = ctx->eps;
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: /*@
476:    NEPNLEIGSGetEPS - Retrieve the eigensolver object (EPS) associated
477:    to the nonlinear eigenvalue solver.

479:    Collective

481:    Input Parameter:
482: .  nep - nonlinear eigenvalue solver

484:    Output Parameter:
485: .  eps - the eigensolver object

487:    Level: advanced

489: .seealso: NEPNLEIGSSetEPS()
490: @*/
491: PetscErrorCode NEPNLEIGSGetEPS(NEP nep,EPS *eps)
492: {
493:   PetscFunctionBegin;
495:   PetscAssertPointer(eps,2);
496:   PetscUseMethod(nep,"NEPNLEIGSGetEPS_C",(NEP,EPS*),(nep,eps));
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }