Actual source code: peprefine.c
slepc-3.19.2 2023-09-05
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 PEP, simple version
12: */
14: #include <slepc/private/pepimpl.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: } PEPSimpNRefctx;
25: typedef struct {
26: Mat M1;
27: Vec M2,M3;
28: PetscScalar M4,m3;
29: } PEP_REFINES_MATSHELL;
31: static PetscErrorCode MatMult_FS(Mat M ,Vec x,Vec y)
32: {
33: PEP_REFINES_MATSHELL *ctx;
34: PetscScalar t;
36: PetscFunctionBegin;
37: PetscCall(MatShellGetContext(M,&ctx));
38: PetscCall(VecDot(x,ctx->M3,&t));
39: t *= ctx->m3/ctx->M4;
40: PetscCall(MatMult(ctx->M1,x,y));
41: PetscCall(VecAXPY(y,-t,ctx->M2));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode PEPSimpleNRefSetUp(PEP pep,PEPSimpNRefctx **ctx_)
46: {
47: PetscInt i,si,j,n0,m0,nloc,*idx1,*idx2,ne;
48: IS is1,is2;
49: PEPSimpNRefctx *ctx;
50: Vec v;
51: PetscMPIInt rank,size;
52: MPI_Comm child;
54: PetscFunctionBegin;
55: PetscCall(PetscCalloc1(1,ctx_));
56: ctx = *ctx_;
57: if (pep->npart==1) {
58: pep->refinesubc = NULL;
59: ctx->scatter_id = NULL;
60: ctx->A = pep->A;
61: } else {
62: PetscCall(PetscSubcommGetChild(pep->refinesubc,&child));
63: PetscCall(PetscMalloc2(pep->nmat,&ctx->A,pep->npart,&ctx->scatter_id));
65: /* Duplicate matrices */
66: for (i=0;i<pep->nmat;i++) PetscCall(MatCreateRedundantMatrix(pep->A[i],0,child,MAT_INITIAL_MATRIX,&ctx->A[i]));
67: PetscCall(MatCreateVecs(ctx->A[0],&ctx->v,NULL));
69: /* Create scatters for sending vectors to each subcommucator */
70: PetscCall(BVGetColumn(pep->V,0,&v));
71: PetscCall(VecGetOwnershipRange(v,&n0,&m0));
72: PetscCall(BVRestoreColumn(pep->V,0,&v));
73: PetscCall(VecGetLocalSize(ctx->v,&nloc));
74: PetscCall(PetscMalloc2(m0-n0,&idx1,m0-n0,&idx2));
75: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),nloc,PETSC_DECIDE,&ctx->vg));
76: for (si=0;si<pep->npart;si++) {
77: j = 0;
78: for (i=n0;i<m0;i++) {
79: idx1[j] = i;
80: idx2[j++] = i+pep->n*si;
81: }
82: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
83: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),(m0-n0),idx2,PETSC_COPY_VALUES,&is2));
84: PetscCall(BVGetColumn(pep->V,0,&v));
85: PetscCall(VecScatterCreate(v,is1,ctx->vg,is2,&ctx->scatter_id[si]));
86: PetscCall(BVRestoreColumn(pep->V,0,&v));
87: PetscCall(ISDestroy(&is1));
88: PetscCall(ISDestroy(&is2));
89: }
90: PetscCall(PetscFree2(idx1,idx2));
91: }
92: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
93: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ctx->A[0]),&rank));
94: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->A[0]),&size));
95: if (size>1) {
96: if (pep->npart==1) PetscCall(BVGetColumn(pep->V,0,&v));
97: else v = ctx->v;
98: PetscCall(VecGetOwnershipRange(v,&n0,&m0));
99: ne = (rank == size-1)?pep->n:0;
100: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->A[0]),ne,PETSC_DECIDE,&ctx->nv));
101: PetscCall(PetscMalloc1(m0-n0,&idx1));
102: for (i=n0;i<m0;i++) idx1[i-n0] = i;
103: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ctx->A[0]),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
104: PetscCall(VecScatterCreate(v,is1,ctx->nv,is1,&ctx->nst));
105: if (pep->npart==1) PetscCall(BVRestoreColumn(pep->V,0,&v));
106: PetscCall(PetscFree(idx1));
107: PetscCall(ISDestroy(&is1));
108: }
109: }
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /*
114: Gather Eigenpair idx from subcommunicator with color sc
115: */
116: static PetscErrorCode PEPSimpleNRefGatherEigenpair(PEP pep,PEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx,PetscInt *fail)
117: {
118: PetscMPIInt nproc,p;
119: MPI_Comm comm=((PetscObject)pep)->comm;
120: Vec v;
121: const PetscScalar *array;
123: PetscFunctionBegin;
124: PetscCallMPI(MPI_Comm_size(comm,&nproc));
125: p = (nproc/pep->npart)*(sc+1)+PetscMin(nproc%pep->npart,sc+1)-1;
126: if (pep->npart>1) {
127: /* Communicate convergence successful */
128: PetscCallMPI(MPI_Bcast(fail,1,MPIU_INT,p,comm));
129: if (!(*fail)) {
130: /* Process 0 of subcommunicator sc broadcasts the eigenvalue */
131: PetscCallMPI(MPI_Bcast(&pep->eigr[idx],1,MPIU_SCALAR,p,comm));
132: /* Gather pep->V[idx] from the subcommuniator sc */
133: PetscCall(BVGetColumn(pep->V,idx,&v));
134: if (pep->refinesubc->color==sc) {
135: PetscCall(VecGetArrayRead(ctx->v,&array));
136: PetscCall(VecPlaceArray(ctx->vg,array));
137: }
138: PetscCall(VecScatterBegin(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
139: PetscCall(VecScatterEnd(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
140: if (pep->refinesubc->color==sc) {
141: PetscCall(VecResetArray(ctx->vg));
142: PetscCall(VecRestoreArrayRead(ctx->v,&array));
143: }
144: PetscCall(BVRestoreColumn(pep->V,idx,&v));
145: }
146: } else {
147: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT && !(*fail)) PetscCallMPI(MPI_Bcast(&pep->eigr[idx],1,MPIU_SCALAR,p,comm));
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode PEPSimpleNRefScatterEigenvector(PEP pep,PEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx)
153: {
154: Vec v;
155: const PetscScalar *array;
157: PetscFunctionBegin;
158: if (pep->npart>1) {
159: PetscCall(BVGetColumn(pep->V,idx,&v));
160: if (pep->refinesubc->color==sc) {
161: PetscCall(VecGetArrayRead(ctx->v,&array));
162: PetscCall(VecPlaceArray(ctx->vg,array));
163: }
164: PetscCall(VecScatterBegin(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
165: PetscCall(VecScatterEnd(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
166: if (pep->refinesubc->color==sc) {
167: PetscCall(VecResetArray(ctx->vg));
168: PetscCall(VecRestoreArrayRead(ctx->v,&array));
169: }
170: PetscCall(BVRestoreColumn(pep->V,idx,&v));
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode PEPEvaluateFunctionDerivatives(PEP pep,PetscScalar alpha,PetscScalar *vals)
176: {
177: PetscInt i,nmat=pep->nmat;
178: PetscScalar a0,a1,a2;
179: PetscReal *a=pep->pbc,*b=a+nmat,*g=b+nmat;
181: PetscFunctionBegin;
182: a0 = 0.0;
183: a1 = 1.0;
184: vals[0] = 0.0;
185: if (nmat>1) vals[1] = 1/a[0];
186: for (i=2;i<nmat;i++) {
187: a2 = ((alpha-b[i-2])*a1-g[i-2]*a0)/a[i-2];
188: vals[i] = (a2+(alpha-b[i-1])*vals[i-1]-g[i-1]*vals[i-2])/a[i-1];
189: a0 = a1; a1 = a2;
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode PEPSimpleNRefSetUpSystem(PEP pep,Mat *A,PEPSimpNRefctx *ctx,PetscInt idx,Mat *Mt,Mat *T,Mat *P,PetscBool ini,Vec t,Vec v)
195: {
196: PetscInt i,nmat=pep->nmat,ml,m0,n0,m1,mg;
197: PetscInt *dnz,*onz,ncols,*cols2=NULL,*nnz;
198: PetscScalar zero=0.0,*coeffs,*coeffs2;
199: PetscMPIInt rank,size;
200: MPI_Comm comm;
201: const PetscInt *cols;
202: const PetscScalar *vals,*array;
203: MatStructure str;
204: PEP_REFINES_MATSHELL *fctx;
205: PEPRefineScheme scheme=pep->scheme;
206: Vec w=ctx->w;
207: Mat M;
209: PetscFunctionBegin;
210: PetscCall(STGetMatStructure(pep->st,&str));
211: PetscCall(PetscMalloc2(nmat,&coeffs,nmat,&coeffs2));
212: switch (scheme) {
213: case PEP_REFINE_SCHEME_SCHUR:
214: if (ini) {
215: PetscCall(PetscCalloc1(1,&fctx));
216: PetscCall(MatGetSize(A[0],&m0,&n0));
217: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A[0]),PETSC_DECIDE,PETSC_DECIDE,m0,n0,fctx,T));
218: PetscCall(MatShellSetOperation(*T,MATOP_MULT,(void(*)(void))MatMult_FS));
219: } else PetscCall(MatShellGetContext(*T,&fctx));
220: M=fctx->M1;
221: break;
222: case PEP_REFINE_SCHEME_MBE:
223: M=*T;
224: break;
225: case PEP_REFINE_SCHEME_EXPLICIT:
226: M=*Mt;
227: break;
228: }
229: if (ini) PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&M));
230: else PetscCall(MatCopy(A[0],M,DIFFERENT_NONZERO_PATTERN));
231: PetscCall(PEPEvaluateBasis(pep,pep->eigr[idx],0,coeffs,NULL));
232: PetscCall(MatScale(M,coeffs[0]));
233: for (i=1;i<nmat;i++) PetscCall(MatAXPY(M,coeffs[i],A[i],(ini)?str:SUBSET_NONZERO_PATTERN));
234: PetscCall(PEPEvaluateFunctionDerivatives(pep,pep->eigr[idx],coeffs2));
235: for (i=0;i<nmat && PetscAbsScalar(coeffs2[i])==0.0;i++);
236: PetscCall(MatMult(A[i],v,w));
237: if (coeffs2[i]!=1.0) PetscCall(VecScale(w,coeffs2[i]));
238: for (i++;i<nmat;i++) {
239: PetscCall(MatMult(A[i],v,t));
240: PetscCall(VecAXPY(w,coeffs2[i],t));
241: }
242: switch (scheme) {
243: case PEP_REFINE_SCHEME_EXPLICIT:
244: comm = PetscObjectComm((PetscObject)A[0]);
245: PetscCallMPI(MPI_Comm_rank(comm,&rank));
246: PetscCallMPI(MPI_Comm_size(comm,&size));
247: PetscCall(MatGetSize(M,&mg,NULL));
248: PetscCall(MatGetOwnershipRange(M,&m0,&m1));
249: if (ini) {
250: PetscCall(MatCreate(comm,T));
251: PetscCall(MatGetLocalSize(M,&ml,NULL));
252: if (rank==size-1) ml++;
253: PetscCall(MatSetSizes(*T,ml,ml,mg+1,mg+1));
254: PetscCall(MatSetFromOptions(*T));
255: PetscCall(MatSetUp(*T));
256: /* Preallocate M */
257: if (size>1) {
258: MatPreallocateBegin(comm,ml,ml,dnz,onz);
259: for (i=m0;i<m1;i++) {
260: PetscCall(MatGetRow(M,i,&ncols,&cols,NULL));
261: PetscCall(MatPreallocateSet(i,ncols,cols,dnz,onz));
262: PetscCall(MatPreallocateSet(i,1,&mg,dnz,onz));
263: PetscCall(MatRestoreRow(M,i,&ncols,&cols,NULL));
264: }
265: if (rank==size-1) {
266: PetscCall(PetscCalloc1(mg+1,&cols2));
267: for (i=0;i<mg+1;i++) cols2[i]=i;
268: PetscCall(MatPreallocateSet(m1,mg+1,cols2,dnz,onz));
269: PetscCall(PetscFree(cols2));
270: }
271: PetscCall(MatMPIAIJSetPreallocation(*T,0,dnz,0,onz));
272: MatPreallocateEnd(dnz,onz);
273: } else {
274: PetscCall(PetscCalloc1(mg+1,&nnz));
275: for (i=0;i<mg;i++) {
276: PetscCall(MatGetRow(M,i,&ncols,NULL,NULL));
277: nnz[i] = ncols+1;
278: PetscCall(MatRestoreRow(M,i,&ncols,NULL,NULL));
279: }
280: nnz[mg] = mg+1;
281: PetscCall(MatSeqAIJSetPreallocation(*T,0,nnz));
282: PetscCall(PetscFree(nnz));
283: }
284: *Mt = M;
285: *P = *T;
286: }
288: /* Set values */
289: PetscCall(VecGetArrayRead(w,&array));
290: for (i=m0;i<m1;i++) {
291: PetscCall(MatGetRow(M,i,&ncols,&cols,&vals));
292: PetscCall(MatSetValues(*T,1,&i,ncols,cols,vals,INSERT_VALUES));
293: PetscCall(MatRestoreRow(M,i,&ncols,&cols,&vals));
294: PetscCall(MatSetValues(*T,1,&i,1,&mg,array+i-m0,INSERT_VALUES));
295: }
296: PetscCall(VecRestoreArrayRead(w,&array));
297: PetscCall(VecConjugate(v));
298: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A[0]),&size));
299: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A[0]),&rank));
300: if (size>1) {
301: if (rank==size-1) {
302: PetscCall(PetscMalloc1(pep->n,&cols2));
303: for (i=0;i<pep->n;i++) cols2[i]=i;
304: }
305: PetscCall(VecScatterBegin(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
306: PetscCall(VecScatterEnd(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
307: PetscCall(VecGetArrayRead(ctx->nv,&array));
308: if (rank==size-1) {
309: PetscCall(MatSetValues(*T,1,&mg,pep->n,cols2,array,INSERT_VALUES));
310: PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
311: }
312: PetscCall(VecRestoreArrayRead(ctx->nv,&array));
313: } else {
314: PetscCall(PetscMalloc1(m1-m0,&cols2));
315: for (i=0;i<m1-m0;i++) cols2[i]=m0+i;
316: PetscCall(VecGetArrayRead(v,&array));
317: PetscCall(MatSetValues(*T,1,&mg,m1-m0,cols2,array,INSERT_VALUES));
318: PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
319: PetscCall(VecRestoreArrayRead(v,&array));
320: }
321: PetscCall(VecConjugate(v));
322: PetscCall(MatAssemblyBegin(*T,MAT_FINAL_ASSEMBLY));
323: PetscCall(MatAssemblyEnd(*T,MAT_FINAL_ASSEMBLY));
324: PetscCall(PetscFree(cols2));
325: break;
326: case PEP_REFINE_SCHEME_SCHUR:
327: fctx->M2 = w;
328: fctx->M3 = v;
329: fctx->m3 = 0.0;
330: for (i=1;i<nmat-1;i++) fctx->m3 += PetscConj(coeffs[i])*coeffs[i];
331: fctx->M4 = 0.0;
332: for (i=1;i<nmat-1;i++) fctx->M4 += PetscConj(coeffs[i])*coeffs2[i];
333: fctx->M1 = M;
334: if (ini) PetscCall(MatDuplicate(M,MAT_COPY_VALUES,P));
335: else PetscCall(MatCopy(M,*P,SAME_NONZERO_PATTERN));
336: if (fctx->M4!=0.0) {
337: PetscCall(VecConjugate(v));
338: PetscCall(VecPointwiseMult(t,v,w));
339: PetscCall(VecConjugate(v));
340: PetscCall(VecScale(t,-fctx->m3/fctx->M4));
341: PetscCall(MatDiagonalSet(*P,t,ADD_VALUES));
342: }
343: break;
344: case PEP_REFINE_SCHEME_MBE:
345: *T = M;
346: *P = M;
347: break;
348: }
349: PetscCall(PetscFree2(coeffs,coeffs2));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: PetscErrorCode PEPNewtonRefinementSimple(PEP pep,PetscInt *maxits,PetscReal tol,PetscInt k)
354: {
355: PetscInt i,n,its,idx=0,*idx_sc,*its_sc,color,*fail_sc;
356: PetscMPIInt rank,size;
357: Mat Mt=NULL,T=NULL,P=NULL;
358: MPI_Comm comm;
359: Vec r,v,dv,rr=NULL,dvv=NULL,t[2];
360: PetscScalar *array2,deig=0.0,tt[2],ttt;
361: const PetscScalar *array;
362: PetscReal norm,error;
363: PetscBool ini=PETSC_TRUE,sc_pend,solved=PETSC_FALSE;
364: PEPSimpNRefctx *ctx;
365: PEP_REFINES_MATSHELL *fctx=NULL;
366: KSPConvergedReason reason;
368: PetscFunctionBegin;
369: PetscCall(PetscLogEventBegin(PEP_Refine,pep,0,0,0));
370: PetscCall(PEPSimpleNRefSetUp(pep,&ctx));
371: its = (maxits)?*maxits:NREF_MAXIT;
372: if (!pep->refineksp) PetscCall(PEPRefineGetKSP(pep,&pep->refineksp));
373: if (pep->npart==1) PetscCall(BVGetColumn(pep->V,0,&v));
374: else v = ctx->v;
375: PetscCall(VecDuplicate(v,&ctx->w));
376: PetscCall(VecDuplicate(v,&r));
377: PetscCall(VecDuplicate(v,&dv));
378: PetscCall(VecDuplicate(v,&t[0]));
379: PetscCall(VecDuplicate(v,&t[1]));
380: if (pep->npart==1) {
381: PetscCall(BVRestoreColumn(pep->V,0,&v));
382: PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
383: } else PetscCall(PetscSubcommGetChild(pep->refinesubc,&comm));
384: PetscCallMPI(MPI_Comm_size(comm,&size));
385: PetscCallMPI(MPI_Comm_rank(comm,&rank));
386: PetscCall(VecGetLocalSize(r,&n));
387: PetscCall(PetscMalloc3(pep->npart,&idx_sc,pep->npart,&its_sc,pep->npart,&fail_sc));
388: for (i=0;i<pep->npart;i++) fail_sc[i] = 0;
389: for (i=0;i<pep->npart;i++) its_sc[i] = 0;
390: color = (pep->npart==1)?0:pep->refinesubc->color;
392: /* Loop performing iterative refinements */
393: while (!solved) {
394: for (i=0;i<pep->npart;i++) {
395: sc_pend = PETSC_TRUE;
396: if (its_sc[i]==0) {
397: idx_sc[i] = idx++;
398: if (idx_sc[i]>=k) {
399: sc_pend = PETSC_FALSE;
400: } else PetscCall(PEPSimpleNRefScatterEigenvector(pep,ctx,i,idx_sc[i]));
401: } else { /* Gather Eigenpair from subcommunicator i */
402: PetscCall(PEPSimpleNRefGatherEigenpair(pep,ctx,i,idx_sc[i],&fail_sc[i]));
403: }
404: while (sc_pend) {
405: if (!fail_sc[i]) PetscCall(PEPComputeError(pep,idx_sc[i],PEP_ERROR_BACKWARD,&error));
406: if (error<=tol || its_sc[i]>=its || fail_sc[i]) {
407: idx_sc[i] = idx++;
408: its_sc[i] = 0;
409: fail_sc[i] = 0;
410: if (idx_sc[i]<k) PetscCall(PEPSimpleNRefScatterEigenvector(pep,ctx,i,idx_sc[i]));
411: } else {
412: sc_pend = PETSC_FALSE;
413: its_sc[i]++;
414: }
415: if (idx_sc[i]>=k) sc_pend = PETSC_FALSE;
416: }
417: }
418: solved = PETSC_TRUE;
419: for (i=0;i<pep->npart&&solved;i++) solved = PetscNot(idx_sc[i]<k);
420: if (idx_sc[color]<k) {
421: #if !defined(PETSC_USE_COMPLEX)
422: PetscCheck(pep->eigi[idx_sc[color]]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Simple Refinement not implemented in real scalars for complex eigenvalues");
423: #endif
424: if (pep->npart==1) PetscCall(BVGetColumn(pep->V,idx_sc[color],&v));
425: else v = ctx->v;
426: PetscCall(PEPSimpleNRefSetUpSystem(pep,ctx->A,ctx,idx_sc[color],&Mt,&T,&P,ini,t[0],v));
427: PetscCall(PEP_KSPSetOperators(pep->refineksp,T,P));
428: if (ini) {
429: PetscCall(KSPSetFromOptions(pep->refineksp));
430: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
431: PetscCall(MatCreateVecs(T,&dvv,NULL));
432: PetscCall(VecDuplicate(dvv,&rr));
433: }
434: ini = PETSC_FALSE;
435: }
437: switch (pep->scheme) {
438: case PEP_REFINE_SCHEME_EXPLICIT:
439: PetscCall(MatMult(Mt,v,r));
440: PetscCall(VecGetArrayRead(r,&array));
441: if (rank==size-1) {
442: PetscCall(VecGetArray(rr,&array2));
443: PetscCall(PetscArraycpy(array2,array,n));
444: array2[n] = 0.0;
445: PetscCall(VecRestoreArray(rr,&array2));
446: } else PetscCall(VecPlaceArray(rr,array));
447: PetscCall(KSPSolve(pep->refineksp,rr,dvv));
448: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
449: if (reason>0) {
450: if (rank != size-1) PetscCall(VecResetArray(rr));
451: PetscCall(VecRestoreArrayRead(r,&array));
452: PetscCall(VecGetArrayRead(dvv,&array));
453: PetscCall(VecPlaceArray(dv,array));
454: PetscCall(VecAXPY(v,-1.0,dv));
455: PetscCall(VecNorm(v,NORM_2,&norm));
456: PetscCall(VecScale(v,1.0/norm));
457: PetscCall(VecResetArray(dv));
458: if (rank==size-1) pep->eigr[idx_sc[color]] -= array[n];
459: PetscCall(VecRestoreArrayRead(dvv,&array));
460: } else fail_sc[color] = 1;
461: break;
462: case PEP_REFINE_SCHEME_MBE:
463: PetscCall(MatMult(T,v,r));
464: /* Mixed block elimination */
465: PetscCall(VecConjugate(v));
466: PetscCall(KSPSolveTranspose(pep->refineksp,v,t[0]));
467: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
468: if (reason>0) {
469: PetscCall(VecConjugate(t[0]));
470: PetscCall(VecDot(ctx->w,t[0],&tt[0]));
471: PetscCall(KSPSolve(pep->refineksp,ctx->w,t[1]));
472: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
473: if (reason>0) {
474: PetscCall(VecDot(t[1],v,&tt[1]));
475: PetscCall(VecDot(r,t[0],&ttt));
476: tt[0] = ttt/tt[0];
477: PetscCall(VecAXPY(r,-tt[0],ctx->w));
478: PetscCall(KSPSolve(pep->refineksp,r,dv));
479: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
480: if (reason>0) {
481: PetscCall(VecDot(dv,v,&ttt));
482: tt[1] = ttt/tt[1];
483: PetscCall(VecAXPY(dv,-tt[1],t[1]));
484: deig = tt[0]+tt[1];
485: }
486: }
487: PetscCall(VecConjugate(v));
488: PetscCall(VecAXPY(v,-1.0,dv));
489: PetscCall(VecNorm(v,NORM_2,&norm));
490: PetscCall(VecScale(v,1.0/norm));
491: pep->eigr[idx_sc[color]] -= deig;
492: fail_sc[color] = 0;
493: } else {
494: PetscCall(VecConjugate(v));
495: fail_sc[color] = 1;
496: }
497: break;
498: case PEP_REFINE_SCHEME_SCHUR:
499: fail_sc[color] = 1;
500: PetscCall(MatShellGetContext(T,&fctx));
501: if (fctx->M4!=0.0) {
502: PetscCall(MatMult(fctx->M1,v,r));
503: PetscCall(KSPSolve(pep->refineksp,r,dv));
504: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
505: if (reason>0) {
506: PetscCall(VecDot(dv,v,&deig));
507: deig *= -fctx->m3/fctx->M4;
508: PetscCall(VecAXPY(v,-1.0,dv));
509: PetscCall(VecNorm(v,NORM_2,&norm));
510: PetscCall(VecScale(v,1.0/norm));
511: pep->eigr[idx_sc[color]] -= deig;
512: fail_sc[color] = 0;
513: }
514: }
515: break;
516: }
517: if (pep->npart==1) PetscCall(BVRestoreColumn(pep->V,idx_sc[color],&v));
518: }
519: }
520: PetscCall(VecDestroy(&t[0]));
521: PetscCall(VecDestroy(&t[1]));
522: PetscCall(VecDestroy(&dv));
523: PetscCall(VecDestroy(&ctx->w));
524: PetscCall(VecDestroy(&r));
525: PetscCall(PetscFree3(idx_sc,its_sc,fail_sc));
526: PetscCall(VecScatterDestroy(&ctx->nst));
527: if (pep->npart>1) {
528: PetscCall(VecDestroy(&ctx->vg));
529: PetscCall(VecDestroy(&ctx->v));
530: for (i=0;i<pep->nmat;i++) PetscCall(MatDestroy(&ctx->A[i]));
531: for (i=0;i<pep->npart;i++) PetscCall(VecScatterDestroy(&ctx->scatter_id[i]));
532: PetscCall(PetscFree2(ctx->A,ctx->scatter_id));
533: }
534: if (fctx && pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
535: PetscCall(MatDestroy(&P));
536: PetscCall(MatDestroy(&fctx->M1));
537: PetscCall(PetscFree(fctx));
538: }
539: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
540: PetscCall(MatDestroy(&Mt));
541: PetscCall(VecDestroy(&dvv));
542: PetscCall(VecDestroy(&rr));
543: PetscCall(VecDestroy(&ctx->nv));
544: }
545: PetscCall(MatDestroy(&T));
546: PetscCall(PetscFree(ctx));
547: PetscCall(PetscLogEventEnd(PEP_Refine,pep,0,0,0));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }