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