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