Actual source code: peprefine.c
slepc-3.21.1 2024-04-26
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 ncols,*cols2=NULL;
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: *Mt = M;
256: *P = *T;
257: }
259: /* Set values */
260: PetscCall(VecGetArrayRead(w,&array));
261: for (i=m0;i<m1;i++) {
262: PetscCall(MatGetRow(M,i,&ncols,&cols,&vals));
263: PetscCall(MatSetValues(*T,1,&i,ncols,cols,vals,INSERT_VALUES));
264: PetscCall(MatRestoreRow(M,i,&ncols,&cols,&vals));
265: PetscCall(MatSetValues(*T,1,&i,1,&mg,array+i-m0,INSERT_VALUES));
266: }
267: PetscCall(VecRestoreArrayRead(w,&array));
268: PetscCall(VecConjugate(v));
269: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A[0]),&size));
270: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A[0]),&rank));
271: if (size>1) {
272: if (rank==size-1) {
273: PetscCall(PetscMalloc1(pep->n,&cols2));
274: for (i=0;i<pep->n;i++) cols2[i]=i;
275: }
276: PetscCall(VecScatterBegin(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
277: PetscCall(VecScatterEnd(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
278: PetscCall(VecGetArrayRead(ctx->nv,&array));
279: if (rank==size-1) {
280: PetscCall(MatSetValues(*T,1,&mg,pep->n,cols2,array,INSERT_VALUES));
281: PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
282: }
283: PetscCall(VecRestoreArrayRead(ctx->nv,&array));
284: } else {
285: PetscCall(PetscMalloc1(m1-m0,&cols2));
286: for (i=0;i<m1-m0;i++) cols2[i]=m0+i;
287: PetscCall(VecGetArrayRead(v,&array));
288: PetscCall(MatSetValues(*T,1,&mg,m1-m0,cols2,array,INSERT_VALUES));
289: PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
290: PetscCall(VecRestoreArrayRead(v,&array));
291: }
292: PetscCall(VecConjugate(v));
293: PetscCall(MatAssemblyBegin(*T,MAT_FINAL_ASSEMBLY));
294: PetscCall(MatAssemblyEnd(*T,MAT_FINAL_ASSEMBLY));
295: PetscCall(PetscFree(cols2));
296: break;
297: case PEP_REFINE_SCHEME_SCHUR:
298: fctx->M2 = w;
299: fctx->M3 = v;
300: fctx->m3 = 0.0;
301: for (i=1;i<nmat-1;i++) fctx->m3 += PetscConj(coeffs[i])*coeffs[i];
302: fctx->M4 = 0.0;
303: for (i=1;i<nmat-1;i++) fctx->M4 += PetscConj(coeffs[i])*coeffs2[i];
304: fctx->M1 = M;
305: if (ini) PetscCall(MatDuplicate(M,MAT_COPY_VALUES,P));
306: else PetscCall(MatCopy(M,*P,SAME_NONZERO_PATTERN));
307: if (fctx->M4!=0.0) {
308: PetscCall(VecConjugate(v));
309: PetscCall(VecPointwiseMult(t,v,w));
310: PetscCall(VecConjugate(v));
311: PetscCall(VecScale(t,-fctx->m3/fctx->M4));
312: PetscCall(MatDiagonalSet(*P,t,ADD_VALUES));
313: }
314: break;
315: case PEP_REFINE_SCHEME_MBE:
316: *T = M;
317: *P = M;
318: break;
319: }
320: PetscCall(PetscFree2(coeffs,coeffs2));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: PetscErrorCode PEPNewtonRefinementSimple(PEP pep,PetscInt *maxits,PetscReal tol,PetscInt k)
325: {
326: PetscInt i,n,its,idx=0,*idx_sc,*its_sc,color,*fail_sc;
327: PetscMPIInt rank,size;
328: Mat Mt=NULL,T=NULL,P=NULL;
329: MPI_Comm comm;
330: Vec r,v,dv,rr=NULL,dvv=NULL,t[2];
331: PetscScalar *array2,deig=0.0,tt[2],ttt;
332: const PetscScalar *array;
333: PetscReal norm,error;
334: PetscBool ini=PETSC_TRUE,sc_pend,solved=PETSC_FALSE;
335: PEPSimpNRefctx *ctx;
336: PEP_REFINES_MATSHELL *fctx=NULL;
337: KSPConvergedReason reason;
339: PetscFunctionBegin;
340: PetscCall(PetscLogEventBegin(PEP_Refine,pep,0,0,0));
341: PetscCall(PEPSimpleNRefSetUp(pep,&ctx));
342: its = (maxits)?*maxits:NREF_MAXIT;
343: if (!pep->refineksp) PetscCall(PEPRefineGetKSP(pep,&pep->refineksp));
344: if (pep->npart==1) PetscCall(BVGetColumn(pep->V,0,&v));
345: else v = ctx->v;
346: PetscCall(VecDuplicate(v,&ctx->w));
347: PetscCall(VecDuplicate(v,&r));
348: PetscCall(VecDuplicate(v,&dv));
349: PetscCall(VecDuplicate(v,&t[0]));
350: PetscCall(VecDuplicate(v,&t[1]));
351: if (pep->npart==1) {
352: PetscCall(BVRestoreColumn(pep->V,0,&v));
353: PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
354: } else PetscCall(PetscSubcommGetChild(pep->refinesubc,&comm));
355: PetscCallMPI(MPI_Comm_size(comm,&size));
356: PetscCallMPI(MPI_Comm_rank(comm,&rank));
357: PetscCall(VecGetLocalSize(r,&n));
358: PetscCall(PetscMalloc3(pep->npart,&idx_sc,pep->npart,&its_sc,pep->npart,&fail_sc));
359: for (i=0;i<pep->npart;i++) fail_sc[i] = 0;
360: for (i=0;i<pep->npart;i++) its_sc[i] = 0;
361: color = (pep->npart==1)?0:pep->refinesubc->color;
363: /* Loop performing iterative refinements */
364: while (!solved) {
365: for (i=0;i<pep->npart;i++) {
366: sc_pend = PETSC_TRUE;
367: if (its_sc[i]==0) {
368: idx_sc[i] = idx++;
369: if (idx_sc[i]>=k) {
370: sc_pend = PETSC_FALSE;
371: } else PetscCall(PEPSimpleNRefScatterEigenvector(pep,ctx,i,idx_sc[i]));
372: } else { /* Gather Eigenpair from subcommunicator i */
373: PetscCall(PEPSimpleNRefGatherEigenpair(pep,ctx,i,idx_sc[i],&fail_sc[i]));
374: }
375: while (sc_pend) {
376: if (!fail_sc[i]) PetscCall(PEPComputeError(pep,idx_sc[i],PEP_ERROR_BACKWARD,&error));
377: if (error<=tol || its_sc[i]>=its || fail_sc[i]) {
378: idx_sc[i] = idx++;
379: its_sc[i] = 0;
380: fail_sc[i] = 0;
381: if (idx_sc[i]<k) PetscCall(PEPSimpleNRefScatterEigenvector(pep,ctx,i,idx_sc[i]));
382: } else {
383: sc_pend = PETSC_FALSE;
384: its_sc[i]++;
385: }
386: if (idx_sc[i]>=k) sc_pend = PETSC_FALSE;
387: }
388: }
389: solved = PETSC_TRUE;
390: for (i=0;i<pep->npart&&solved;i++) solved = PetscNot(idx_sc[i]<k);
391: if (idx_sc[color]<k) {
392: #if !defined(PETSC_USE_COMPLEX)
393: PetscCheck(pep->eigi[idx_sc[color]]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Simple Refinement not implemented in real scalars for complex eigenvalues");
394: #endif
395: if (pep->npart==1) PetscCall(BVGetColumn(pep->V,idx_sc[color],&v));
396: else v = ctx->v;
397: PetscCall(PEPSimpleNRefSetUpSystem(pep,ctx->A,ctx,idx_sc[color],&Mt,&T,&P,ini,t[0],v));
398: PetscCall(PEP_KSPSetOperators(pep->refineksp,T,P));
399: if (ini) {
400: PetscCall(KSPSetFromOptions(pep->refineksp));
401: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
402: PetscCall(MatCreateVecs(T,&dvv,NULL));
403: PetscCall(VecDuplicate(dvv,&rr));
404: }
405: ini = PETSC_FALSE;
406: }
408: switch (pep->scheme) {
409: case PEP_REFINE_SCHEME_EXPLICIT:
410: PetscCall(MatMult(Mt,v,r));
411: PetscCall(VecGetArrayRead(r,&array));
412: if (rank==size-1) {
413: PetscCall(VecGetArray(rr,&array2));
414: PetscCall(PetscArraycpy(array2,array,n));
415: array2[n] = 0.0;
416: PetscCall(VecRestoreArray(rr,&array2));
417: } else PetscCall(VecPlaceArray(rr,array));
418: PetscCall(KSPSolve(pep->refineksp,rr,dvv));
419: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
420: if (reason>0) {
421: if (rank != size-1) PetscCall(VecResetArray(rr));
422: PetscCall(VecRestoreArrayRead(r,&array));
423: PetscCall(VecGetArrayRead(dvv,&array));
424: PetscCall(VecPlaceArray(dv,array));
425: PetscCall(VecAXPY(v,-1.0,dv));
426: PetscCall(VecNorm(v,NORM_2,&norm));
427: PetscCall(VecScale(v,1.0/norm));
428: PetscCall(VecResetArray(dv));
429: if (rank==size-1) pep->eigr[idx_sc[color]] -= array[n];
430: PetscCall(VecRestoreArrayRead(dvv,&array));
431: } else fail_sc[color] = 1;
432: break;
433: case PEP_REFINE_SCHEME_MBE:
434: PetscCall(MatMult(T,v,r));
435: /* Mixed block elimination */
436: PetscCall(VecConjugate(v));
437: PetscCall(KSPSolveTranspose(pep->refineksp,v,t[0]));
438: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
439: if (reason>0) {
440: PetscCall(VecConjugate(t[0]));
441: PetscCall(VecDot(ctx->w,t[0],&tt[0]));
442: PetscCall(KSPSolve(pep->refineksp,ctx->w,t[1]));
443: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
444: if (reason>0) {
445: PetscCall(VecDot(t[1],v,&tt[1]));
446: PetscCall(VecDot(r,t[0],&ttt));
447: tt[0] = ttt/tt[0];
448: PetscCall(VecAXPY(r,-tt[0],ctx->w));
449: PetscCall(KSPSolve(pep->refineksp,r,dv));
450: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
451: if (reason>0) {
452: PetscCall(VecDot(dv,v,&ttt));
453: tt[1] = ttt/tt[1];
454: PetscCall(VecAXPY(dv,-tt[1],t[1]));
455: deig = tt[0]+tt[1];
456: }
457: }
458: PetscCall(VecConjugate(v));
459: PetscCall(VecAXPY(v,-1.0,dv));
460: PetscCall(VecNorm(v,NORM_2,&norm));
461: PetscCall(VecScale(v,1.0/norm));
462: pep->eigr[idx_sc[color]] -= deig;
463: fail_sc[color] = 0;
464: } else {
465: PetscCall(VecConjugate(v));
466: fail_sc[color] = 1;
467: }
468: break;
469: case PEP_REFINE_SCHEME_SCHUR:
470: fail_sc[color] = 1;
471: PetscCall(MatShellGetContext(T,&fctx));
472: if (fctx->M4!=0.0) {
473: PetscCall(MatMult(fctx->M1,v,r));
474: PetscCall(KSPSolve(pep->refineksp,r,dv));
475: PetscCall(KSPGetConvergedReason(pep->refineksp,&reason));
476: if (reason>0) {
477: PetscCall(VecDot(dv,v,&deig));
478: deig *= -fctx->m3/fctx->M4;
479: PetscCall(VecAXPY(v,-1.0,dv));
480: PetscCall(VecNorm(v,NORM_2,&norm));
481: PetscCall(VecScale(v,1.0/norm));
482: pep->eigr[idx_sc[color]] -= deig;
483: fail_sc[color] = 0;
484: }
485: }
486: break;
487: }
488: if (pep->npart==1) PetscCall(BVRestoreColumn(pep->V,idx_sc[color],&v));
489: }
490: }
491: PetscCall(VecDestroy(&t[0]));
492: PetscCall(VecDestroy(&t[1]));
493: PetscCall(VecDestroy(&dv));
494: PetscCall(VecDestroy(&ctx->w));
495: PetscCall(VecDestroy(&r));
496: PetscCall(PetscFree3(idx_sc,its_sc,fail_sc));
497: PetscCall(VecScatterDestroy(&ctx->nst));
498: if (pep->npart>1) {
499: PetscCall(VecDestroy(&ctx->vg));
500: PetscCall(VecDestroy(&ctx->v));
501: for (i=0;i<pep->nmat;i++) PetscCall(MatDestroy(&ctx->A[i]));
502: for (i=0;i<pep->npart;i++) PetscCall(VecScatterDestroy(&ctx->scatter_id[i]));
503: PetscCall(PetscFree2(ctx->A,ctx->scatter_id));
504: }
505: if (fctx && pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
506: PetscCall(MatDestroy(&P));
507: PetscCall(MatDestroy(&fctx->M1));
508: PetscCall(PetscFree(fctx));
509: }
510: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
511: PetscCall(MatDestroy(&Mt));
512: PetscCall(VecDestroy(&dvv));
513: PetscCall(VecDestroy(&rr));
514: PetscCall(VecDestroy(&ctx->nv));
515: }
516: PetscCall(MatDestroy(&T));
517: PetscCall(PetscFree(ctx));
518: PetscCall(PetscLogEventEnd(PEP_Refine,pep,0,0,0));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }