Actual source code: neprefine.c
slepc-3.20.2 2024-03-15
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: }