Actual source code: slp-twosided.c

slepc-3.16.0 2021-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    Two-sided variant of the NEPSLP solver.
 12: */

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

 17: typedef struct _n_nep_def_ctx *NEP_NEDEF_CTX;

 19: struct _n_nep_def_ctx {
 20:   PetscInt    n;
 21:   PetscBool   ref;
 22:   PetscScalar *eig;
 23:   BV          V,W;
 24: };

 26: typedef struct {   /* context for two-sided solver */
 27:   Mat         Ft;
 28:   Mat         Jt;
 29:   Vec         w;
 30: } NEP_SLPTS_MATSHELL;

 32: typedef struct {   /* context for non-equivalence deflation */
 33:   NEP_NEDEF_CTX defctx;
 34:   Mat           F;
 35:   Mat           J;
 36:   KSP           ksp;
 37:   PetscBool     isJ;
 38:   PetscScalar   lambda;
 39:   Vec           w[2];
 40: } NEP_NEDEF_MATSHELL;

 42: static PetscErrorCode MatMult_SLPTS_Right(Mat M,Vec x,Vec y)
 43: {
 44:   PetscErrorCode     ierr;
 45:   NEP_SLPTS_MATSHELL *ctx;

 48:   MatShellGetContext(M,&ctx);
 49:   MatMult(ctx->Jt,x,ctx->w);
 50:   MatSolve(ctx->Ft,ctx->w,y);
 51:   return(0);
 52: }

 54: static PetscErrorCode MatMult_SLPTS_Left(Mat M,Vec x,Vec y)
 55: {
 56:   PetscErrorCode     ierr;
 57:   NEP_SLPTS_MATSHELL *ctx;

 60:   MatShellGetContext(M,&ctx);
 61:   MatMultTranspose(ctx->Jt,x,ctx->w);
 62:   MatSolveTranspose(ctx->Ft,ctx->w,y);
 63:   return(0);
 64: }

 66: static PetscErrorCode MatDestroy_SLPTS(Mat M)
 67: {
 68:   PetscErrorCode     ierr;
 69:   NEP_SLPTS_MATSHELL *ctx;

 72:   MatShellGetContext(M,&ctx);
 73:   VecDestroy(&ctx->w);
 74:   PetscFree(ctx);
 75:   return(0);
 76: }

 78: #if defined(PETSC_HAVE_CUDA)
 79: static PetscErrorCode MatCreateVecs_SLPTS(Mat M,Vec *left,Vec *right)
 80: {
 81:   PetscErrorCode     ierr;
 82:   NEP_SLPTS_MATSHELL *ctx;

 85:   MatShellGetContext(M,&ctx);
 86:   if (right) {
 87:     VecDuplicate(ctx->w,right);
 88:   }
 89:   if (left) {
 90:     VecDuplicate(ctx->w,left);
 91:   }
 92:   return(0);
 93: }
 94: #endif

 96: static PetscErrorCode NEPSLPSetUpEPSMat(NEP nep,Mat F,Mat J,PetscBool left,Mat *M)
 97: {
 98:   PetscErrorCode     ierr;
 99:   Mat                Mshell;
100:   PetscInt           nloc,mloc;
101:   NEP_SLPTS_MATSHELL *shellctx;

104:   /* Create mat shell */
105:   PetscNew(&shellctx);
106:   shellctx->Ft = F;
107:   shellctx->Jt = J;
108:   MatGetLocalSize(nep->function,&mloc,&nloc);
109:   MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell);
110:   if (left) {
111:     MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Left);
112:   } else {
113:     MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Right);
114:   }
115:   MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_SLPTS);
116: #if defined(PETSC_HAVE_CUDA)
117:   MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_SLPTS);
118: #endif
119:   *M = Mshell;
120:   MatCreateVecs(nep->function,&shellctx->w,NULL);
121:   return(0);
122: }

124: /* Functions for deflation */
125: static PetscErrorCode NEPDeflationNEDestroy(NEP_NEDEF_CTX defctx)
126: {

130:   if (!defctx) return(0);
131:   PetscFree(defctx->eig);
132:   PetscFree(defctx);
133:   return(0);
134: }

136: static PetscErrorCode NEPDeflationNECreate(NEP nep,BV V,BV W,PetscInt sz,NEP_NEDEF_CTX *defctx)
137: {
139:   NEP_NEDEF_CTX  op;

142:   PetscNew(&op);
143:   *defctx = op;
144:   op->n   = 0;
145:   op->ref = PETSC_FALSE;
146:   PetscCalloc1(sz,&op->eig);
147:   PetscObjectStateIncrease((PetscObject)V);
148:   PetscObjectStateIncrease((PetscObject)W);
149:   op->V = V;
150:   op->W = W;
151:   return(0);
152: }

154: static PetscErrorCode NEPDeflationNEComputeFunction(NEP nep,Mat M,PetscScalar lambda)
155: {
156:   PetscErrorCode     ierr;
157:   NEP_NEDEF_MATSHELL *matctx;

160:   MatShellGetContext(M,&matctx);
161:   if (lambda==matctx->lambda) return(0);
162:   NEPComputeFunction(nep,lambda,matctx->F,matctx->F);
163:   if (matctx->isJ) {NEPComputeJacobian(nep,lambda,matctx->J);}
164:   if (matctx->ksp) {KSPSetOperators(matctx->ksp,matctx->F,matctx->F);}
165:   matctx->lambda = lambda;
166:   return(0);
167: }

169: static PetscErrorCode MatMult_NEPDeflationNE(Mat M,Vec x,Vec r)
170: {
171:   Vec                t,tt;
172:   PetscScalar        *h,*alpha,lambda,*eig;
173:   PetscInt           i,k;
174:   PetscErrorCode     ierr;
175:   NEP_NEDEF_MATSHELL *matctx;

178:   MatShellGetContext(M,&matctx);
179:   if (matctx->defctx->n && !matctx->defctx->ref) {
180:     k = matctx->defctx->n;
181:     lambda = matctx->lambda;
182:     eig = matctx->defctx->eig;
183:     t = matctx->w[0];
184:     VecCopy(x,t);
185:     PetscMalloc2(k,&h,k,&alpha);
186:     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
187:     BVDotVec(matctx->defctx->V,t,h);
188:     for (i=0;i<k;i++) h[i] *= alpha[i];
189:     BVMultVec(matctx->defctx->W,-1.0,1.0,t,h);
190:     MatMult(matctx->isJ?matctx->J:matctx->F,t,r);
191:     if (matctx->isJ) {
192:       for (i=0;i<k;i++) h[i] *= (1.0/((lambda-eig[i])*(lambda-eig[i])))/alpha[i];
193:       tt = matctx->w[1];
194:       BVMultVec(matctx->defctx->W,-1.0,0.0,tt,h);
195:       MatMult(matctx->F,tt,t);
196:       VecAXPY(r,1.0,t);
197:     }
198:     PetscFree2(h,alpha);
199:   } else {
200:     MatMult(matctx->isJ?matctx->J:matctx->F,x,r);
201:   }
202:   return(0);
203: }

205: static PetscErrorCode MatMultTranspose_NEPDeflationNE(Mat M,Vec x,Vec r)
206: {
207:   Vec                t,tt;
208:   PetscScalar        *h,*alphaC,lambda,*eig;
209:   PetscInt           i,k;
210:   PetscErrorCode     ierr;
211:   NEP_NEDEF_MATSHELL *matctx;

214:   MatShellGetContext(M,&matctx);
215:   t    = matctx->w[0];
216:   VecCopy(x,t);
217:   if (matctx->defctx->n && !matctx->defctx->ref) {
218:     VecConjugate(t);
219:     k = matctx->defctx->n;
220:     lambda = matctx->lambda;
221:     eig = matctx->defctx->eig;
222:     PetscMalloc2(k,&h,k,&alphaC);
223:     for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
224:     BVDotVec(matctx->defctx->W,t,h);
225:     for (i=0;i<k;i++) h[i] *= alphaC[i];
226:     BVMultVec(matctx->defctx->V,-1.0,1.0,t,h);
227:     VecConjugate(t);
228:     MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r);
229:     if (matctx->isJ) {
230:       for (i=0;i<k;i++) h[i] *= PetscConj(1.0/((lambda-eig[i])*(lambda-eig[i])))/alphaC[i];
231:       tt = matctx->w[1];
232:       BVMultVec(matctx->defctx->V,-1.0,0.0,tt,h);
233:       VecConjugate(tt);
234:       MatMultTranspose(matctx->F,tt,t);
235:       VecAXPY(r,1.0,t);
236:     }
237:     PetscFree2(h,alphaC);
238:   } else {
239:     MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r);
240:   }
241:   return(0);
242: }

244: static PetscErrorCode MatSolve_NEPDeflationNE(Mat M,Vec b,Vec x)
245: {
246:   PetscErrorCode     ierr;
247:   PetscScalar        *h,*alpha,lambda,*eig;
248:   PetscInt           i,k;
249:   NEP_NEDEF_MATSHELL *matctx;

252:   MatShellGetContext(M,&matctx);
253:   if (!matctx->ksp) {
254:     VecCopy(b,x);
255:     return(0);
256:   }
257:   KSPSolve(matctx->ksp,b,x);
258:   if (matctx->defctx->n && !matctx->defctx->ref) {
259:     k = matctx->defctx->n;
260:     lambda = matctx->lambda;
261:     eig = matctx->defctx->eig;
262:     PetscMalloc2(k,&h,k,&alpha);
263:     BVDotVec(matctx->defctx->V,x,h);
264:     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
265:     for (i=0;i<k;i++) h[i] *= alpha[i]/(1.0-alpha[i]);
266:     BVMultVec(matctx->defctx->W,1.0,1.0,x,h);
267:     PetscFree2(h,alpha);
268:   }
269:   return(0);
270: }

272: static PetscErrorCode MatSolveTranspose_NEPDeflationNE(Mat M,Vec b,Vec x)
273: {
274:   PetscErrorCode     ierr;
275:   PetscScalar        *h,*alphaC,lambda,*eig;
276:   PetscInt           i,k;
277:   NEP_NEDEF_MATSHELL *matctx;

280:   MatShellGetContext(M,&matctx);
281:   if (!matctx->ksp) {
282:     VecCopy(b,x);
283:     return(0);
284:   }
285:   KSPSolveTranspose(matctx->ksp,b,x);
286:   if (matctx->defctx->n && !matctx->defctx->ref) {
287:     VecConjugate(x);
288:     k = matctx->defctx->n;
289:     lambda = matctx->lambda;
290:     eig = matctx->defctx->eig;
291:     PetscMalloc2(k,&h,k,&alphaC);
292:     BVDotVec(matctx->defctx->W,x,h);
293:     for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
294:     for (i=0;i<k;i++) h[i] *= alphaC[i]/(1.0-alphaC[i]);
295:     BVMultVec(matctx->defctx->V,1.0,1.0,x,h);
296:     PetscFree2(h,alphaC);
297:     VecConjugate(x);
298:   }
299:   return(0);
300: }

302: static PetscErrorCode MatDestroy_NEPDeflationNE(Mat M)
303: {
304:   PetscErrorCode     ierr;
305:   NEP_NEDEF_MATSHELL *matctx;

308:   MatShellGetContext(M,&matctx);
309:   VecDestroy(&matctx->w[0]);
310:   VecDestroy(&matctx->w[1]);
311:   PetscFree(matctx);
312:   return(0);
313: }

315: static PetscErrorCode MatCreateVecs_NEPDeflationNE(Mat M,Vec *right,Vec *left)
316: {
317:   PetscErrorCode     ierr;
318:   NEP_NEDEF_MATSHELL *matctx;

321:   MatShellGetContext(M,&matctx);
322:   MatCreateVecs(matctx->F,right,left);
323:   return(0);
324: }

326: static PetscErrorCode NEPDeflationNEFunctionCreate(NEP_NEDEF_CTX defctx,NEP nep,Mat F,Mat J,KSP ksp,PetscBool isJ,Mat *Mshell)
327: {
328:   NEP_NEDEF_MATSHELL *matctx;
329:   PetscErrorCode     ierr;
330:   PetscInt           nloc,mloc;

333:   /* Create mat shell */
334:   PetscNew(&matctx);
335:   MatGetLocalSize(nep->function,&mloc,&nloc);
336:   MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Mshell);
337:   matctx->F   = F;
338:   matctx->J   = J;
339:   matctx->isJ = isJ;
340:   matctx->ksp = ksp;
341:   matctx->defctx = defctx;
342:   matctx->lambda = PETSC_MAX_REAL;
343:   MatCreateVecs(F,&matctx->w[0],NULL);
344:   VecDuplicate(matctx->w[0],&matctx->w[1]);
345:   MatShellSetOperation(*Mshell,MATOP_MULT,(void(*)(void))MatMult_NEPDeflationNE);
346:   MatShellSetOperation(*Mshell,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_NEPDeflationNE);
347:   MatShellSetOperation(*Mshell,MATOP_SOLVE,(void(*)(void))MatSolve_NEPDeflationNE);
348:   MatShellSetOperation(*Mshell,MATOP_SOLVE_TRANSPOSE,(void(*)(void))MatSolveTranspose_NEPDeflationNE);
349:   MatShellSetOperation(*Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_NEPDeflationNE);
350:   MatShellSetOperation(*Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_NEPDeflationNE);
351:   return(0);
352: }

354: static PetscErrorCode NEPDeflationNERecoverEigenvectors(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
355: {
357:   PetscScalar    *h,*alpha,*eig;
358:   PetscInt       i,k;

361:   if (w) { VecConjugate(w); }
362:   if (defctx->n && !defctx->ref) {
363:     eig = defctx->eig;
364:     k = defctx->n;
365:     PetscMalloc2(k,&h,k,&alpha);
366:     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
367:     BVDotVec(defctx->V,u,h);
368:     for (i=0;i<k;i++) h[i] *= alpha[i];
369:     BVMultVec(defctx->W,-1.0,1.0,u,h);
370:     VecNormalize(u,NULL);
371:     if (w) {
372:       BVDotVec(defctx->W,w,h);
373:       for (i=0;i<k;i++) alpha[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
374:       for (i=0;i<k;i++) h[i] *= alpha[i];
375:       BVMultVec(defctx->V,-1.0,1.0,w,h);
376:       VecNormalize(w,NULL);
377:     }
378:     PetscFree2(h,alpha);
379:   }
380:   return(0);
381: }

383: static PetscErrorCode NEPDeflationNELocking(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
384: {
386:   PetscInt       n;

389:   n = defctx->n++;
390:   defctx->eig[n] = lambda;
391:   BVInsertVec(defctx->V,n,u);
392:   BVInsertVec(defctx->W,n,w);
393:   BVSetActiveColumns(defctx->V,0,defctx->n);
394:   BVSetActiveColumns(defctx->W,0,defctx->n);
395:   BVBiorthonormalizeColumn(defctx->V,defctx->W,n,NULL);
396:   return(0);
397: }

399: static PetscErrorCode NEPDeflationNESetRefine(NEP_NEDEF_CTX defctx,PetscBool ref)
400: {
402:   defctx->ref = ref;
403:   return(0);
404: }

406: PetscErrorCode NEPSolve_SLP_Twosided(NEP nep)
407: {
409:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
410:   Mat            mF,mJ,M,Mt;
411:   Vec            u,r,t,w;
412:   BV             X,Y;
413:   PetscScalar    sigma,lambda,mu,im=0.0,mu2,im2;
414:   PetscReal      resnorm,resl;
415:   PetscInt       nconv,nconv2,i;
416:   PetscBool      skip=PETSC_FALSE,lock=PETSC_FALSE;
417:   NEP_NEDEF_CTX  defctx=NULL;    /* Extended operator for deflation */

420:   /* get initial approximation of eigenvalue and eigenvector */
421:   NEPGetDefaultShift(nep,&sigma);
422:   if (!nep->nini) {
423:     BVSetRandomColumn(nep->V,0);
424:   }
425:   BVSetRandomColumn(nep->W,0);
426:   lambda = sigma;
427:   if (!ctx->ksp) { NEPSLPGetKSP(nep,&ctx->ksp); }
428:   BVDuplicate(nep->V,&X);
429:   BVDuplicate(nep->W,&Y);
430:   NEPDeflationNECreate(nep,X,Y,nep->nev,&defctx);
431:   BVGetColumn(nep->V,0,&t);
432:   VecDuplicate(t,&u);
433:   VecDuplicate(t,&w);
434:   BVRestoreColumn(nep->V,0,&t);
435:   BVCopyVec(nep->V,0,u);
436:   BVCopyVec(nep->W,0,w);
437:   VecDuplicate(u,&r);
438:   NEPDeflationNEFunctionCreate(defctx,nep,nep->function,NULL,ctx->ksp,PETSC_FALSE,&mF);
439:   NEPDeflationNEFunctionCreate(defctx,nep,nep->function,nep->jacobian,NULL,PETSC_TRUE,&mJ);
440:   NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_FALSE,&M);
441:   NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_TRUE,&Mt);
442:   EPSSetOperators(ctx->eps,M,NULL);
443:   MatDestroy(&M);
444:   EPSSetOperators(ctx->epsts,Mt,NULL);
445:   MatDestroy(&Mt);

447:   /* Restart loop */
448:   while (nep->reason == NEP_CONVERGED_ITERATING) {
449:     nep->its++;

451:     /* form residual,  r = T(lambda)*u (used in convergence test only) */
452:     NEPDeflationNEComputeFunction(nep,mF,lambda);
453:     MatMultTranspose(mF,w,r);
454:     VecNorm(r,NORM_2,&resl);
455:     MatMult(mF,u,r);

457:     /* convergence test */
458:     VecNorm(r,NORM_2,&resnorm);
459:     resnorm = PetscMax(resnorm,resl);
460:     (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
461:     nep->eigr[nep->nconv] = lambda;
462:     if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
463:       if (nep->errest[nep->nconv]<=ctx->deftol && !defctx->ref && nep->nconv) {
464:         NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda);
465:         VecConjugate(w);
466:         NEPDeflationNESetRefine(defctx,PETSC_TRUE);
467:         MatMultTranspose(mF,w,r);
468:         VecNorm(r,NORM_2,&resl);
469:         MatMult(mF,u,r);
470:         VecNorm(r,NORM_2,&resnorm);
471:         resnorm = PetscMax(resnorm,resl);
472:         (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
473:         if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
474:       } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
475:     }
476:     if (lock) {
477:       lock = PETSC_FALSE;
478:       skip = PETSC_TRUE;
479:       NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda);
480:       NEPDeflationNELocking(defctx,u,w,lambda);
481:       NEPDeflationNESetRefine(defctx,PETSC_FALSE);
482:       BVInsertVec(nep->V,nep->nconv,u);
483:       BVInsertVec(nep->W,nep->nconv,w);
484:       VecConjugate(w);
485:       nep->nconv = nep->nconv + 1;
486:     }
487:     (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
488:     if (!skip || nep->reason>0) {
489:       NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
490:     }

492:     if (nep->reason == NEP_CONVERGED_ITERATING) {
493:       if (!skip) {
494:         /* evaluate T(lambda) and T'(lambda) */
495:         NEPDeflationNEComputeFunction(nep,mF,lambda);
496:         NEPDeflationNEComputeFunction(nep,mJ,lambda);
497:         EPSSetInitialSpace(ctx->eps,1,&u);
498:         EPSSetInitialSpace(ctx->epsts,1,&w);

500:         /* compute new eigenvalue correction mu and eigenvector approximation u */
501:         EPSSolve(ctx->eps);
502:         EPSSolve(ctx->epsts);
503:         EPSGetConverged(ctx->eps,&nconv);
504:         EPSGetConverged(ctx->epsts,&nconv2);
505:         if (!nconv||!nconv2) {
506:           PetscInfo1(nep,"iter=%D, inner iteration failed, stopping solve\n",nep->its);
507:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
508:           break;
509:         }
510:         EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL);
511:         for (i=0;i<nconv2;i++) {
512:           EPSGetEigenpair(ctx->epsts,i,&mu2,&im2,w,NULL);
513:           if (SlepcAbsEigenvalue(mu-mu2,im-im2)/SlepcAbsEigenvalue(mu,im)<nep->tol*1000) break;
514:         }
515:         if (i==nconv2) {
516:           PetscInfo1(nep,"iter=%D, inner iteration failed, stopping solve\n",nep->its);
517:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
518:           break;
519:         }

521:         mu = 1.0/mu;
522:         if (PetscAbsScalar(im)>PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
523:       } else {
524:         nep->its--;  /* do not count this as a full iteration */
525:         /* use second eigenpair computed in previous iteration */
526:         EPSGetConverged(ctx->eps,&nconv);
527:         if (nconv>=2 && nconv2>=2) {
528:           EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL);
529:           EPSGetEigenpair(ctx->epsts,1,&mu2,&im2,w,NULL);
530:           mu = 1.0/mu;
531:         } else {
532:           BVSetRandomColumn(nep->V,nep->nconv);
533:           BVSetRandomColumn(nep->W,nep->nconv);
534:           BVCopyVec(nep->V,nep->nconv,u);
535:           BVCopyVec(nep->W,nep->nconv,w);
536:           mu = lambda-sigma;
537:         }
538:         skip = PETSC_FALSE;
539:       }
540:       /* correct eigenvalue */
541:       lambda = lambda - mu;
542:     }
543:   }
544:   VecDestroy(&u);
545:   VecDestroy(&w);
546:   VecDestroy(&r);
547:   MatDestroy(&mF);
548:   MatDestroy(&mJ);
549:   BVDestroy(&X);
550:   BVDestroy(&Y);
551:   NEPDeflationNEDestroy(defctx);
552:   return(0);
553: }