Actual source code: neprefine.c

slepc-3.19.2 2023-09-05
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:    Newton refinement for NEP, simple version
 12: */

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

 17: #define NREF_MAXIT 10

 19: typedef struct {
 20:   VecScatter    *scatter_id,nst;
 21:   Mat           *A;
 22:   Vec           nv,vg,v,w;
 23:   FN            *fn;
 24: } NEPSimpNRefctx;

 26: typedef struct {
 27:   Mat          M1;
 28:   Vec          M2,M3;
 29:   PetscScalar  M4,m3;
 30: } NEP_REFINE_MATSHELL;

 32: static PetscErrorCode MatMult_FS(Mat M ,Vec x,Vec y)
 33: {
 34:   NEP_REFINE_MATSHELL *ctx;
 35:   PetscScalar         t;

 37:   PetscFunctionBegin;
 38:   PetscCall(MatShellGetContext(M,&ctx));
 39:   PetscCall(VecDot(x,ctx->M3,&t));
 40:   t *= ctx->m3/ctx->M4;
 41:   PetscCall(MatMult(ctx->M1,x,y));
 42:   PetscCall(VecAXPY(y,-t,ctx->M2));
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: static PetscErrorCode NEPSimpleNRefSetUp(NEP nep,NEPSimpNRefctx **ctx_)
 47: {
 48:   PetscInt       i,si,j,n0,m0,nloc,*idx1,*idx2,ne;
 49:   IS             is1,is2;
 50:   NEPSimpNRefctx *ctx;
 51:   Vec            v;
 52:   PetscMPIInt    rank,size;
 53:   MPI_Comm       child;

 55:   PetscFunctionBegin;
 56:   PetscCall(PetscCalloc1(1,ctx_));
 57:   ctx = *ctx_;
 58:   if (nep->npart==1) {
 59:     ctx->A  = nep->A;
 60:     ctx->fn = nep->f;
 61:     nep->refinesubc = NULL;
 62:     ctx->scatter_id = NULL;
 63:   } else {
 64:     PetscCall(PetscSubcommGetChild(nep->refinesubc,&child));
 65:     PetscCall(PetscMalloc2(nep->nt,&ctx->A,nep->npart,&ctx->scatter_id));

 67:     /* Duplicate matrices */
 68:     for (i=0;i<nep->nt;i++) PetscCall(MatCreateRedundantMatrix(nep->A[i],0,child,MAT_INITIAL_MATRIX,&ctx->A[i]));
 69:     PetscCall(MatCreateVecs(ctx->A[0],&ctx->v,NULL));

 71:     /* Duplicate FNs */
 72:     PetscCall(PetscMalloc1(nep->nt,&ctx->fn));
 73:     for (i=0;i<nep->nt;i++) PetscCall(FNDuplicate(nep->f[i],child,&ctx->fn[i]));

 75:     /* Create scatters for sending vectors to each subcommucator */
 76:     PetscCall(BVGetColumn(nep->V,0,&v));
 77:     PetscCall(VecGetOwnershipRange(v,&n0,&m0));
 78:     PetscCall(BVRestoreColumn(nep->V,0,&v));
 79:     PetscCall(VecGetLocalSize(ctx->v,&nloc));
 80:     PetscCall(PetscMalloc2(m0-n0,&idx1,m0-n0,&idx2));
 81:     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)nep),nloc,PETSC_DECIDE,&ctx->vg));
 82:     for (si=0;si<nep->npart;si++) {
 83:       j = 0;
 84:       for (i=n0;i<m0;i++) {
 85:         idx1[j]   = i;
 86:         idx2[j++] = i+nep->n*si;
 87:       }
 88:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)nep),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
 89:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)nep),(m0-n0),idx2,PETSC_COPY_VALUES,&is2));
 90:       PetscCall(BVGetColumn(nep->V,0,&v));
 91:       PetscCall(VecScatterCreate(v,is1,ctx->vg,is2,&ctx->scatter_id[si]));
 92:       PetscCall(BVRestoreColumn(nep->V,0,&v));
 93:       PetscCall(ISDestroy(&is1));
 94:       PetscCall(ISDestroy(&is2));
 95:     }
 96:     PetscCall(PetscFree2(idx1,idx2));
 97:   }
 98:   if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
 99:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ctx->A[0]),&rank));
100:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->A[0]),&size));
101:     if (size>1) {
102:       if (nep->npart==1) PetscCall(BVGetColumn(nep->V,0,&v));
103:       else v = ctx->v;
104:       PetscCall(VecGetOwnershipRange(v,&n0,&m0));
105:       ne = (rank == size-1)?nep->n:0;
106:       PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->A[0]),ne,PETSC_DECIDE,&ctx->nv));
107:       PetscCall(PetscMalloc1(m0-n0,&idx1));
108:       for (i=n0;i<m0;i++) {
109:         idx1[i-n0] = i;
110:       }
111:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ctx->A[0]),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
112:       PetscCall(VecScatterCreate(v,is1,ctx->nv,is1,&ctx->nst));
113:       if (nep->npart==1) PetscCall(BVRestoreColumn(nep->V,0,&v));
114:       PetscCall(PetscFree(idx1));
115:       PetscCall(ISDestroy(&is1));
116:     }
117:   }  PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /*
121:   Gather Eigenpair idx from subcommunicator with color sc
122: */
123: static PetscErrorCode NEPSimpleNRefGatherEigenpair(NEP nep,NEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx,PetscInt *fail)
124: {
125:   PetscMPIInt    nproc,p;
126:   MPI_Comm       comm=((PetscObject)nep)->comm;
127:   Vec            v;
128:   PetscScalar    *array;

130:   PetscFunctionBegin;
131:   PetscCallMPI(MPI_Comm_size(comm,&nproc));
132:   p = (nproc/nep->npart)*(sc+1)+PetscMin(sc+1,nproc%nep->npart)-1;
133:   if (nep->npart>1) {
134:     /* Communicate convergence successful */
135:     PetscCallMPI(MPI_Bcast(fail,1,MPIU_INT,p,comm));
136:     if (!(*fail)) {
137:       /* Process 0 of subcommunicator sc broadcasts the eigenvalue */
138:       PetscCallMPI(MPI_Bcast(&nep->eigr[idx],1,MPIU_SCALAR,p,comm));
139:       /* Gather nep->V[idx] from the subcommuniator sc */
140:       PetscCall(BVGetColumn(nep->V,idx,&v));
141:       if (nep->refinesubc->color==sc) {
142:         PetscCall(VecGetArray(ctx->v,&array));
143:         PetscCall(VecPlaceArray(ctx->vg,array));
144:       }
145:       PetscCall(VecScatterBegin(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
146:       PetscCall(VecScatterEnd(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
147:       if (nep->refinesubc->color==sc) {
148:         PetscCall(VecResetArray(ctx->vg));
149:         PetscCall(VecRestoreArray(ctx->v,&array));
150:       }
151:       PetscCall(BVRestoreColumn(nep->V,idx,&v));
152:     }
153:   } else {
154:     if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT && !(*fail)) PetscCallMPI(MPI_Bcast(&nep->eigr[idx],1,MPIU_SCALAR,p,comm));
155:   }
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode NEPSimpleNRefScatterEigenvector(NEP nep,NEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx)
160: {
161:   Vec            v;
162:   PetscScalar    *array;

164:   PetscFunctionBegin;
165:   if (nep->npart>1) {
166:     PetscCall(BVGetColumn(nep->V,idx,&v));
167:     if (nep->refinesubc->color==sc) {
168:       PetscCall(VecGetArray(ctx->v,&array));
169:       PetscCall(VecPlaceArray(ctx->vg,array));
170:     }
171:     PetscCall(VecScatterBegin(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
172:     PetscCall(VecScatterEnd(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
173:     if (nep->refinesubc->color==sc) {
174:       PetscCall(VecResetArray(ctx->vg));
175:       PetscCall(VecRestoreArray(ctx->v,&array));
176:     }
177:     PetscCall(BVRestoreColumn(nep->V,idx,&v));
178:   }
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

182: static PetscErrorCode NEPSimpleNRefSetUpSystem(NEP nep,NEPSimpNRefctx *ctx,Mat *A,PetscInt idx,Mat *Mt,Mat *T,Mat *P,PetscBool ini,Vec t,Vec v)
183: {
184:   PetscInt            i,st,ml,m0,n0,m1,mg;
185:   PetscInt            *dnz,*onz,ncols,*cols2=NULL,*nnz,nt=nep->nt;
186:   PetscScalar         zero=0.0,*coeffs,*coeffs2;
187:   PetscMPIInt         rank,size;
188:   MPI_Comm            comm;
189:   const PetscInt      *cols;
190:   const PetscScalar   *vals,*array;
191:   NEP_REFINE_MATSHELL *fctx;
192:   Vec                 w=ctx->w;
193:   Mat                 M;

195:   PetscFunctionBegin;
196:   PetscCall(PetscMalloc2(nt,&coeffs,nt,&coeffs2));
197:   switch (nep->scheme) {
198:   case NEP_REFINE_SCHEME_SCHUR:
199:     if (ini) {
200:       PetscCall(PetscCalloc1(1,&fctx));
201:       PetscCall(MatGetSize(A[0],&m0,&n0));
202:       PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A[0]),PETSC_DECIDE,PETSC_DECIDE,m0,n0,fctx,T));
203:       PetscCall(MatShellSetOperation(*T,MATOP_MULT,(void(*)(void))MatMult_FS));
204:     } else PetscCall(MatShellGetContext(*T,&fctx));
205:     M=fctx->M1;
206:     break;
207:   case NEP_REFINE_SCHEME_MBE:
208:     M=*T;
209:     break;
210:   case NEP_REFINE_SCHEME_EXPLICIT:
211:     M=*Mt;
212:     break;
213:   }
214:   if (ini) PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&M));
215:   else PetscCall(MatCopy(A[0],M,DIFFERENT_NONZERO_PATTERN));
216:   for (i=0;i<nt;i++) PetscCall(FNEvaluateFunction(ctx->fn[i],nep->eigr[idx],coeffs+i));
217:   if (coeffs[0]!=1.0) PetscCall(MatScale(M,coeffs[0]));
218:   for (i=1;i<nt;i++) PetscCall(MatAXPY(M,coeffs[i],A[i],(ini)?nep->mstr:SUBSET_NONZERO_PATTERN));
219:   for (i=0;i<nt;i++) PetscCall(FNEvaluateDerivative(ctx->fn[i],nep->eigr[idx],coeffs2+i));
220:   st = 0;
221:   for (i=0;i<nt && PetscAbsScalar(coeffs2[i])==0.0;i++) st++;
222:   PetscCall(MatMult(A[st],v,w));
223:   if (coeffs2[st]!=1.0) PetscCall(VecScale(w,coeffs2[st]));
224:   for (i=st+1;i<nt;i++) {
225:     PetscCall(MatMult(A[i],v,t));
226:     PetscCall(VecAXPY(w,coeffs2[i],t));
227:   }

229:   switch (nep->scheme) {
230:   case NEP_REFINE_SCHEME_EXPLICIT:
231:     comm = PetscObjectComm((PetscObject)A[0]);
232:     PetscCallMPI(MPI_Comm_rank(comm,&rank));
233:     PetscCallMPI(MPI_Comm_size(comm,&size));
234:     PetscCall(MatGetSize(M,&mg,NULL));
235:     PetscCall(MatGetOwnershipRange(M,&m0,&m1));
236:     if (ini) {
237:       PetscCall(MatCreate(comm,T));
238:       PetscCall(MatGetLocalSize(M,&ml,NULL));
239:       if (rank==size-1) ml++;
240:       PetscCall(MatSetSizes(*T,ml,ml,mg+1,mg+1));
241:       PetscCall(MatSetFromOptions(*T));
242:       PetscCall(MatSetUp(*T));
243:       /* Preallocate M */
244:       if (size>1) {
245:         MatPreallocateBegin(comm,ml,ml,dnz,onz);
246:         for (i=m0;i<m1;i++) {
247:           PetscCall(MatGetRow(M,i,&ncols,&cols,NULL));
248:           PetscCall(MatPreallocateSet(i,ncols,cols,dnz,onz));
249:           PetscCall(MatPreallocateSet(i,1,&mg,dnz,onz));
250:           PetscCall(MatRestoreRow(M,i,&ncols,&cols,NULL));
251:         }
252:         if (rank==size-1) {
253:           PetscCall(PetscCalloc1(mg+1,&cols2));
254:           for (i=0;i<mg+1;i++) cols2[i]=i;
255:           PetscCall(MatPreallocateSet(m1,mg+1,cols2,dnz,onz));
256:           PetscCall(PetscFree(cols2));
257:         }
258:         PetscCall(MatMPIAIJSetPreallocation(*T,0,dnz,0,onz));
259:         MatPreallocateEnd(dnz,onz);
260:       } else {
261:         PetscCall(PetscCalloc1(mg+1,&nnz));
262:         for (i=0;i<mg;i++) {
263:           PetscCall(MatGetRow(M,i,&ncols,NULL,NULL));
264:           nnz[i] = ncols+1;
265:           PetscCall(MatRestoreRow(M,i,&ncols,NULL,NULL));
266:         }
267:         nnz[mg] = mg+1;
268:         PetscCall(MatSeqAIJSetPreallocation(*T,0,nnz));
269:         PetscCall(PetscFree(nnz));
270:       }
271:       *Mt = M;
272:       *P  = *T;
273:     }

275:     /* Set values */
276:     PetscCall(VecGetArrayRead(w,&array));
277:     for (i=m0;i<m1;i++) {
278:       PetscCall(MatGetRow(M,i,&ncols,&cols,&vals));
279:       PetscCall(MatSetValues(*T,1,&i,ncols,cols,vals,INSERT_VALUES));
280:       PetscCall(MatRestoreRow(M,i,&ncols,&cols,&vals));
281:       PetscCall(MatSetValues(*T,1,&i,1,&mg,array+i-m0,INSERT_VALUES));
282:     }
283:     PetscCall(VecRestoreArrayRead(w,&array));
284:     PetscCall(VecConjugate(v));
285:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A[0]),&size));
286:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A[0]),&rank));
287:     if (size>1) {
288:       if (rank==size-1) {
289:         PetscCall(PetscMalloc1(nep->n,&cols2));
290:         for (i=0;i<nep->n;i++) cols2[i]=i;
291:       }
292:       PetscCall(VecScatterBegin(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
293:       PetscCall(VecScatterEnd(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
294:       PetscCall(VecGetArrayRead(ctx->nv,&array));
295:       if (rank==size-1) {
296:         PetscCall(MatSetValues(*T,1,&mg,nep->n,cols2,array,INSERT_VALUES));
297:         PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
298:       }
299:       PetscCall(VecRestoreArrayRead(ctx->nv,&array));
300:     } else {
301:       PetscCall(PetscMalloc1(m1-m0,&cols2));
302:       for (i=0;i<m1-m0;i++) cols2[i]=m0+i;
303:       PetscCall(VecGetArrayRead(v,&array));
304:       PetscCall(MatSetValues(*T,1,&mg,m1-m0,cols2,array,INSERT_VALUES));
305:       PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
306:       PetscCall(VecRestoreArrayRead(v,&array));
307:     }
308:     PetscCall(VecConjugate(v));
309:     PetscCall(MatAssemblyBegin(*T,MAT_FINAL_ASSEMBLY));
310:     PetscCall(MatAssemblyEnd(*T,MAT_FINAL_ASSEMBLY));
311:     PetscCall(PetscFree(cols2));
312:     break;
313:   case NEP_REFINE_SCHEME_SCHUR:
314:     fctx->M2 = ctx->w;
315:     fctx->M3 = v;
316:     fctx->m3 = 1.0+PetscConj(nep->eigr[idx])*nep->eigr[idx];
317:     fctx->M4 = PetscConj(nep->eigr[idx]);
318:     fctx->M1 = M;
319:     if (ini) PetscCall(MatDuplicate(M,MAT_COPY_VALUES,P));
320:     else PetscCall(MatCopy(M,*P,SAME_NONZERO_PATTERN));
321:     if (fctx->M4!=0.0) {
322:       PetscCall(VecConjugate(v));
323:       PetscCall(VecPointwiseMult(t,v,w));
324:       PetscCall(VecConjugate(v));
325:       PetscCall(VecScale(t,-fctx->m3/fctx->M4));
326:       PetscCall(MatDiagonalSet(*P,t,ADD_VALUES));
327:     }
328:     break;
329:   case NEP_REFINE_SCHEME_MBE:
330:     *T = M;
331:     *P = M;
332:     break;
333:   }
334:   PetscCall(PetscFree2(coeffs,coeffs2));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: PetscErrorCode NEPNewtonRefinementSimple(NEP nep,PetscInt *maxits,PetscReal tol,PetscInt k)
339: {
340:   PetscInt            i,n,its,idx=0,*idx_sc,*its_sc,color,*fail_sc;
341:   PetscMPIInt         rank,size;
342:   Mat                 Mt=NULL,T=NULL,P=NULL;
343:   MPI_Comm            comm;
344:   Vec                 r,v,dv,rr=NULL,dvv=NULL,t[2];
345:   const PetscScalar   *array;
346:   PetscScalar         *array2,deig=0.0,tt[2],ttt;
347:   PetscReal           norm,error;
348:   PetscBool           ini=PETSC_TRUE,sc_pend,solved=PETSC_FALSE;
349:   NEPSimpNRefctx      *ctx;
350:   NEP_REFINE_MATSHELL *fctx=NULL;
351:   KSPConvergedReason  reason;

353:   PetscFunctionBegin;
354:   PetscCall(PetscLogEventBegin(NEP_Refine,nep,0,0,0));
355:   PetscCall(NEPSimpleNRefSetUp(nep,&ctx));
356:   its = (maxits)?*maxits:NREF_MAXIT;
357:   if (!nep->refineksp) PetscCall(NEPRefineGetKSP(nep,&nep->refineksp));
358:   if (nep->npart==1) PetscCall(BVGetColumn(nep->V,0,&v));
359:   else v = ctx->v;
360:   PetscCall(VecDuplicate(v,&ctx->w));
361:   PetscCall(VecDuplicate(v,&r));
362:   PetscCall(VecDuplicate(v,&dv));
363:   PetscCall(VecDuplicate(v,&t[0]));
364:   PetscCall(VecDuplicate(v,&t[1]));
365:   if (nep->npart==1) {
366:     PetscCall(BVRestoreColumn(nep->V,0,&v));
367:     PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
368:   } else PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
369:   PetscCallMPI(MPI_Comm_size(comm,&size));
370:   PetscCallMPI(MPI_Comm_rank(comm,&rank));
371:   PetscCall(VecGetLocalSize(r,&n));
372:   PetscCall(PetscMalloc3(nep->npart,&idx_sc,nep->npart,&its_sc,nep->npart,&fail_sc));
373:   for (i=0;i<nep->npart;i++) fail_sc[i] = 0;
374:   for (i=0;i<nep->npart;i++) its_sc[i] = 0;
375:   color = (nep->npart==1)?0:nep->refinesubc->color;

377:   /* Loop performing iterative refinements */
378:   while (!solved) {
379:     for (i=0;i<nep->npart;i++) {
380:       sc_pend = PETSC_TRUE;
381:       if (its_sc[i]==0) {
382:         idx_sc[i] = idx++;
383:         if (idx_sc[i]>=k) {
384:           sc_pend = PETSC_FALSE;
385:         } else PetscCall(NEPSimpleNRefScatterEigenvector(nep,ctx,i,idx_sc[i]));
386:       }  else { /* Gather Eigenpair from subcommunicator i */
387:         PetscCall(NEPSimpleNRefGatherEigenpair(nep,ctx,i,idx_sc[i],&fail_sc[i]));
388:       }
389:       while (sc_pend) {
390:         if (!fail_sc[i]) PetscCall(NEPComputeError(nep,idx_sc[i],NEP_ERROR_RELATIVE,&error));
391:         if (error<=tol || its_sc[i]>=its || fail_sc[i]) {
392:           idx_sc[i] = idx++;
393:           its_sc[i] = 0;
394:           fail_sc[i] = 0;
395:           if (idx_sc[i]<k) PetscCall(NEPSimpleNRefScatterEigenvector(nep,ctx,i,idx_sc[i]));
396:         } else {
397:           sc_pend = PETSC_FALSE;
398:           its_sc[i]++;
399:         }
400:         if (idx_sc[i]>=k) sc_pend = PETSC_FALSE;
401:       }
402:     }
403:     solved = PETSC_TRUE;
404:     for (i=0;i<nep->npart&&solved;i++) solved = PetscNot(idx_sc[i]<k);
405:     if (idx_sc[color]<k) {
406: #if !defined(PETSC_USE_COMPLEX)
407:       PetscCheck(nep->eigi[idx_sc[color]]==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Simple Refinement not implemented in real scalar for complex eigenvalues");
408: #endif
409:       if (nep->npart==1) PetscCall(BVGetColumn(nep->V,idx_sc[color],&v));
410:       else v = ctx->v;
411:       PetscCall(NEPSimpleNRefSetUpSystem(nep,ctx,ctx->A,idx_sc[color],&Mt,&T,&P,ini,t[0],v));
412:       PetscCall(NEP_KSPSetOperators(nep->refineksp,T,P));
413:       if (ini) {
414:         PetscCall(KSPSetFromOptions(nep->refineksp));
415:         if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
416:           PetscCall(MatCreateVecs(T,&dvv,NULL));
417:           PetscCall(VecDuplicate(dvv,&rr));
418:         }
419:         ini = PETSC_FALSE;
420:       }
421:       switch (nep->scheme) {
422:       case NEP_REFINE_SCHEME_EXPLICIT:
423:         PetscCall(MatMult(Mt,v,r));
424:         PetscCall(VecGetArrayRead(r,&array));
425:         if (rank==size-1) {
426:           PetscCall(VecGetArray(rr,&array2));
427:           PetscCall(PetscArraycpy(array2,array,n));
428:           array2[n] = 0.0;
429:           PetscCall(VecRestoreArray(rr,&array2));
430:         } else PetscCall(VecPlaceArray(rr,array));
431:         PetscCall(KSPSolve(nep->refineksp,rr,dvv));
432:         PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
433:         if (reason>0) {
434:           if (rank != size-1) PetscCall(VecResetArray(rr));
435:           PetscCall(VecRestoreArrayRead(r,&array));
436:           PetscCall(VecGetArrayRead(dvv,&array));
437:           PetscCall(VecPlaceArray(dv,array));
438:           PetscCall(VecAXPY(v,-1.0,dv));
439:           PetscCall(VecNorm(v,NORM_2,&norm));
440:           PetscCall(VecScale(v,1.0/norm));
441:           PetscCall(VecResetArray(dv));
442:           if (rank==size-1) nep->eigr[idx_sc[color]] -= array[n];
443:           PetscCall(VecRestoreArrayRead(dvv,&array));
444:         } else fail_sc[color] = 1;
445:         break;
446:       case NEP_REFINE_SCHEME_MBE:
447:         PetscCall(MatMult(T,v,r));
448:         /* Mixed block elimination */
449:         PetscCall(VecConjugate(v));
450:         PetscCall(KSPSolveTranspose(nep->refineksp,v,t[0]));
451:         PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
452:         if (reason>0) {
453:           PetscCall(VecConjugate(t[0]));
454:           PetscCall(VecDot(ctx->w,t[0],&tt[0]));
455:           PetscCall(KSPSolve(nep->refineksp,ctx->w,t[1]));
456:           PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
457:           if (reason>0) {
458:             PetscCall(VecDot(t[1],v,&tt[1]));
459:             PetscCall(VecDot(r,t[0],&ttt));
460:             tt[0] = ttt/tt[0];
461:             PetscCall(VecAXPY(r,-tt[0],ctx->w));
462:             PetscCall(KSPSolve(nep->refineksp,r,dv));
463:             PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
464:             if (reason>0) {
465:               PetscCall(VecDot(dv,v,&ttt));
466:               tt[1] = ttt/tt[1];
467:               PetscCall(VecAXPY(dv,-tt[1],t[1]));
468:               deig = tt[0]+tt[1];
469:             }
470:           }
471:           PetscCall(VecConjugate(v));
472:           PetscCall(VecAXPY(v,-1.0,dv));
473:           PetscCall(VecNorm(v,NORM_2,&norm));
474:           PetscCall(VecScale(v,1.0/norm));
475:           nep->eigr[idx_sc[color]] -= deig;
476:           fail_sc[color] = 0;
477:         } else {
478:           PetscCall(VecConjugate(v));
479:           fail_sc[color] = 1;
480:         }
481:         break;
482:       case NEP_REFINE_SCHEME_SCHUR:
483:         fail_sc[color] = 1;
484:         PetscCall(MatShellGetContext(T,&fctx));
485:         if (fctx->M4!=0.0) {
486:           PetscCall(MatMult(fctx->M1,v,r));
487:           PetscCall(KSPSolve(nep->refineksp,r,dv));
488:           PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
489:           if (reason>0) {
490:             PetscCall(VecDot(dv,v,&deig));
491:             deig *= -fctx->m3/fctx->M4;
492:             PetscCall(VecAXPY(v,-1.0,dv));
493:             PetscCall(VecNorm(v,NORM_2,&norm));
494:             PetscCall(VecScale(v,1.0/norm));
495:             nep->eigr[idx_sc[color]] -= deig;
496:             fail_sc[color] = 0;
497:           }
498:         }
499:         break;
500:       }
501:       if (nep->npart==1) PetscCall(BVRestoreColumn(nep->V,idx_sc[color],&v));
502:     }
503:   }
504:   PetscCall(VecDestroy(&t[0]));
505:   PetscCall(VecDestroy(&t[1]));
506:   PetscCall(VecDestroy(&dv));
507:   PetscCall(VecDestroy(&ctx->w));
508:   PetscCall(VecDestroy(&r));
509:   PetscCall(PetscFree3(idx_sc,its_sc,fail_sc));
510:   PetscCall(VecScatterDestroy(&ctx->nst));
511:   if (nep->npart>1) {
512:     PetscCall(VecDestroy(&ctx->vg));
513:     PetscCall(VecDestroy(&ctx->v));
514:     for (i=0;i<nep->nt;i++) PetscCall(MatDestroy(&ctx->A[i]));
515:     for (i=0;i<nep->npart;i++) PetscCall(VecScatterDestroy(&ctx->scatter_id[i]));
516:     PetscCall(PetscFree2(ctx->A,ctx->scatter_id));
517:   }
518:   if (fctx && nep->scheme==NEP_REFINE_SCHEME_SCHUR) {
519:     PetscCall(MatDestroy(&P));
520:     PetscCall(MatDestroy(&fctx->M1));
521:     PetscCall(PetscFree(fctx));
522:   }
523:   if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
524:     PetscCall(MatDestroy(&Mt));
525:     PetscCall(VecDestroy(&dvv));
526:     PetscCall(VecDestroy(&rr));
527:     PetscCall(VecDestroy(&ctx->nv));
528:     if (nep->npart>1) {
529:       for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&ctx->fn[i]));
530:       PetscCall(PetscFree(ctx->fn));
531:     }
532:   }
533:   if (nep->scheme==NEP_REFINE_SCHEME_MBE) {
534:     if (nep->npart>1) {
535:       for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&ctx->fn[i]));
536:       PetscCall(PetscFree(ctx->fn));
537:     }
538:   }
539:   PetscCall(MatDestroy(&T));
540:   PetscCall(PetscFree(ctx));
541:   PetscCall(PetscLogEventEnd(NEP_Refine,nep,0,0,0));
542:   PetscFunctionReturn(PETSC_SUCCESS);
543: }