Actual source code: nrefine.c
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 polynomial eigenproblems.
13: References:
15: [1] T. Betcke and D. Kressner, "Perturbation, extraction and refinement
16: of invariant pairs for matrix polynomials", Linear Algebra Appl.
17: 435(3):514-536, 2011.
19: [2] C. Campos and J.E. Roman, "Parallel iterative refinement in
20: polynomial eigenvalue problems", Numer. Linear Algebra Appl. 23(4):
21: 730-745, 2016.
22: */
24: #include <slepc/private/pepimpl.h>
25: #include <slepcblaslapack.h>
27: typedef struct {
28: Mat *A,M1;
29: BV V,M2,M3,W;
30: PetscInt k,nmat;
31: PetscScalar *fih,*work,*M4;
32: PetscBLASInt *pM4;
33: PetscBool compM1;
34: Vec t;
35: } PEP_REFINE_MATSHELL;
37: typedef struct {
38: Mat E[2],M1;
39: Vec tN,ttN,t1,vseq;
40: VecScatter scatterctx;
41: PetscBool compM1;
42: PetscInt *map0,*map1,*idxg,*idxp;
43: PetscSubcomm subc;
44: VecScatter scatter_sub;
45: VecScatter *scatter_id,*scatterp_id;
46: Mat *A;
47: BV V,W,M2,M3,Wt;
48: PetscScalar *M4,*w,*wt,*d,*dt;
49: Vec t,tg,Rv,Vi,tp,tpg;
50: PetscInt idx,*cols;
51: } PEP_REFINE_EXPLICIT;
53: static PetscErrorCode MatMult_FS(Mat M ,Vec x,Vec y)
54: {
55: PEP_REFINE_MATSHELL *ctx;
56: PetscInt k,i;
57: PetscScalar *c;
58: PetscBLASInt k_,one=1;
60: PetscFunctionBegin;
61: PetscCall(MatShellGetContext(M,&ctx));
62: PetscCall(VecCopy(x,ctx->t));
63: k = ctx->k;
64: c = ctx->work;
65: PetscCall(PetscBLASIntCast(k,&k_));
66: PetscCall(MatMult(ctx->M1,x,y));
67: PetscCall(VecConjugate(ctx->t));
68: PetscCall(BVDotVec(ctx->M3,ctx->t,c));
69: for (i=0;i<k;i++) c[i] = PetscConj(c[i]);
70: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
71: PetscCallLAPACKInfo("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,c,&k_,&info));
72: PetscCall(PetscFPTrapPop());
73: PetscCall(BVMultVec(ctx->M2,-1.0,1.0,y,c));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*
78: Evaluates the first d elements of the polynomial basis
79: on a given matrix H which is considered to be triangular
80: */
81: static PetscErrorCode PEPEvaluateBasisforMatrix(PEP pep,PetscInt nm,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH)
82: {
83: PetscInt i,j,ldfh=nm*k,off,nmat=pep->nmat;
84: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat,t;
85: PetscScalar corr=0.0,alpha,beta;
86: PetscBLASInt k_,ldh_,ldfh_;
88: PetscFunctionBegin;
89: PetscCall(PetscBLASIntCast(ldh,&ldh_));
90: PetscCall(PetscBLASIntCast(k,&k_));
91: PetscCall(PetscBLASIntCast(ldfh,&ldfh_));
92: PetscCall(PetscArrayzero(fH,nm*k*k));
93: if (nm>0) for (j=0;j<k;j++) fH[j+j*ldfh] = 1.0;
94: if (nm>1) {
95: t = b[0]/a[0];
96: off = k;
97: for (j=0;j<k;j++) {
98: for (i=0;i<k;i++) fH[off+i+j*ldfh] = H[i+j*ldh]/a[0];
99: fH[j+j*ldfh] -= t;
100: }
101: }
102: for (i=2;i<nm;i++) {
103: off = i*k;
104: if (i==2) {
105: for (j=0;j<k;j++) {
106: fH[off+j+j*ldfh] = 1.0;
107: H[j+j*ldh] -= b[1];
108: }
109: } else {
110: for (j=0;j<k;j++) {
111: PetscCall(PetscArraycpy(fH+off+j*ldfh,fH+(i-2)*k+j*ldfh,k));
112: H[j+j*ldh] += corr-b[i-1];
113: }
114: }
115: corr = b[i-1];
116: beta = -g[i-1]/a[i-1];
117: alpha = 1/a[i-1];
118: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&alpha,H,&ldh_,fH+(i-1)*k,&ldfh_,&beta,fH+off,&ldfh_));
119: }
120: for (j=0;j<k;j++) H[j+j*ldh] += corr;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode NRefSysSetup_shell(PEP pep,PetscInt k,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,PEP_REFINE_MATSHELL *ctx)
125: {
126: PetscScalar *DHii,*T12,*Tr,*Ts,*array,s,ss,sone=1.0,zero=0.0,*M4=ctx->M4,t,*v,*T;
127: const PetscScalar *m3,*m2;
128: PetscInt i,d,j,nmat=pep->nmat,lda=nmat*k,deg=nmat-1,nloc,ld2,ld3;
129: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
130: PetscBLASInt k_,lda_,lds_,nloc_,ld2_,one=1;
131: Mat *A=ctx->A,Mk,M1=ctx->M1,P;
132: BV V=ctx->V,M2=ctx->M2,M3=ctx->M3,W=ctx->W;
133: MatStructure str;
134: Vec vc;
136: PetscFunctionBegin;
137: PetscCall(STGetMatStructure(pep->st,&str));
138: PetscCall(PetscMalloc3(nmat*k*k,&T12,k*k,&Tr,PetscMax(k*k,nmat),&Ts));
139: DHii = T12;
140: PetscCall(PetscArrayzero(DHii,k*k*nmat));
141: for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
142: for (d=2;d<nmat;d++) {
143: for (j=0;j<k;j++) {
144: for (i=0;i<k;i++) {
145: DHii[d*k+i+j*lda] = ((h-b[d-1])*DHii[(d-1)*k+i+j*lda]+fH[(d-1)*k+i+j*lda]-g[d-1]*DHii[(d-2)*k+i+j*lda])/a[d-1];
146: }
147: }
148: }
149: /* T11 */
150: if (!ctx->compM1) {
151: PetscCall(MatCopy(A[0],M1,DIFFERENT_NONZERO_PATTERN));
152: PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
153: for (j=1;j<nmat;j++) PetscCall(MatAXPY(M1,Ts[j],A[j],str));
154: }
156: /* T22 */
157: PetscCall(PetscBLASIntCast(lds,&lds_));
158: PetscCall(PetscBLASIntCast(k,&k_));
159: PetscCall(PetscBLASIntCast(lda,&lda_));
160: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
161: for (i=1;i<deg;i++) {
162: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
163: s = (i==1)?0.0:1.0;
164: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,M4,&k_));
165: }
166: for (i=0;i<k;i++) for (j=0;j<i;j++) { t=M4[i+j*k];M4[i+j*k]=M4[j+i*k];M4[j+i*k]=t; }
168: /* T12 */
169: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&Mk));
170: for (i=1;i<nmat;i++) {
171: PetscCall(MatDenseGetArrayWrite(Mk,&array));
172: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,array,&k_));
173: PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
174: PetscCall(BVSetActiveColumns(W,0,k));
175: PetscCall(BVMult(W,1.0,0.0,V,Mk));
176: if (i==1) PetscCall(BVMatMult(W,A[i],M2));
177: else {
178: PetscCall(BVMatMult(W,A[i],M3)); /* using M3 as work space */
179: PetscCall(BVMult(M2,1.0,1.0,M3,NULL));
180: }
181: }
183: /* T21 */
184: PetscCall(MatDenseGetArrayWrite(Mk,&array));
185: for (i=1;i<deg;i++) {
186: s = (i==1)?0.0:1.0;
187: ss = PetscConj(fh[i]);
188: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,array,&k_));
189: }
190: PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
191: PetscCall(BVSetActiveColumns(M3,0,k));
192: PetscCall(BVMult(M3,1.0,0.0,V,Mk));
193: for (i=0;i<k;i++) {
194: PetscCall(BVGetColumn(M3,i,&vc));
195: PetscCall(VecConjugate(vc));
196: PetscCall(BVRestoreColumn(M3,i,&vc));
197: }
198: PetscCall(MatDestroy(&Mk));
199: PetscCall(PetscFree3(T12,Tr,Ts));
201: PetscCall(VecGetLocalSize(ctx->t,&nloc));
202: PetscCall(PetscBLASIntCast(nloc,&nloc_));
203: PetscCall(PetscMalloc1(nloc*k,&T));
204: PetscCall(KSPGetOperators(pep->refineksp,NULL,&P));
205: if (!ctx->compM1) PetscCall(MatCopy(ctx->M1,P,SAME_NONZERO_PATTERN));
206: PetscCall(BVGetArrayRead(ctx->M2,&m2));
207: PetscCall(BVGetLeadingDimension(ctx->M2,&ld2));
208: PetscCall(PetscBLASIntCast(ld2,&ld2_));
209: PetscCall(BVGetArrayRead(ctx->M3,&m3));
210: PetscCall(BVGetLeadingDimension(ctx->M3,&ld3));
211: PetscCall(VecGetArray(ctx->t,&v));
212: for (i=0;i<nloc;i++) for (j=0;j<k;j++) T[j+i*k] = m3[i+j*ld3];
213: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
214: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&k_,&nloc_,ctx->M4,&k_,ctx->pM4,T,&k_,&info));
215: PetscCall(PetscFPTrapPop());
216: for (i=0;i<nloc;i++) v[i] = BLASdot_(&k_,m2+i,&ld2_,T+i*k,&one);
217: PetscCall(VecRestoreArray(ctx->t,&v));
218: PetscCall(BVRestoreArrayRead(ctx->M2,&m2));
219: PetscCall(BVRestoreArrayRead(ctx->M3,&m3));
220: PetscCall(MatDiagonalSet(P,ctx->t,ADD_VALUES));
221: PetscCall(PetscFree(T));
222: PetscCall(KSPSetUp(pep->refineksp));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode NRefSysSolve_shell(KSP ksp,PetscInt nmat,Vec Rv,PetscScalar *Rh,PetscInt k,Vec dVi,PetscScalar *dHi)
227: {
228: PetscScalar *t0;
229: PetscBLASInt k_,one=1,lda_;
230: PetscInt i,lda=nmat*k;
231: Mat M;
232: PEP_REFINE_MATSHELL *ctx;
234: PetscFunctionBegin;
235: PetscCall(KSPGetOperators(ksp,&M,NULL));
236: PetscCall(MatShellGetContext(M,&ctx));
237: PetscCall(PetscCalloc1(k,&t0));
238: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
239: PetscCall(PetscBLASIntCast(lda,&lda_));
240: PetscCall(PetscBLASIntCast(k,&k_));
241: for (i=0;i<k;i++) t0[i] = Rh[i];
242: PetscCallLAPACKInfo("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,t0,&k_,&info));
243: PetscCall(BVMultVec(ctx->M2,-1.0,1.0,Rv,t0));
244: PetscCall(KSPSolve(ksp,Rv,dVi));
245: PetscCall(VecConjugate(dVi));
246: PetscCall(BVDotVec(ctx->M3,dVi,dHi));
247: PetscCall(VecConjugate(dVi));
248: for (i=0;i<k;i++) dHi[i] = Rh[i]-PetscConj(dHi[i]);
249: PetscCallLAPACKInfo("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,dHi,&k_,&info));
250: PetscCall(PetscFPTrapPop());
251: PetscCall(PetscFree(t0));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*
256: Computes the residual P(H,V*S)*e_j for the polynomial
257: */
258: static PetscErrorCode NRefRightSide(PetscInt nmat,PetscReal *pcf,Mat *A,PetscInt k,BV V,PetscScalar *S,PetscInt lds,PetscInt j,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *DfH,PetscScalar *dH,BV dV,PetscScalar *dVS,PetscInt rds,Vec Rv,PetscScalar *Rh,BV W,Vec t)
259: {
260: PetscScalar *DS0,*DS1,*F,beta=0.0,sone=1.0,none=-1.0,tt=0.0,*h,zero=0.0,*Z,*c0;
261: PetscReal *a=pcf,*b=pcf+nmat,*g=b+nmat;
262: PetscInt i,ii,jj,lda;
263: PetscBLASInt lda_,k_,ldh_,lds_,nmat_,k2_,krds_,j_,one=1;
264: Mat M0;
265: Vec w;
267: PetscFunctionBegin;
268: PetscCall(PetscMalloc4(k*nmat,&h,k*k,&DS0,k*k,&DS1,k*k,&Z));
269: lda = k*nmat;
270: PetscCall(PetscBLASIntCast(k,&k_));
271: PetscCall(PetscBLASIntCast(lds,&lds_));
272: PetscCall(PetscBLASIntCast(lda,&lda_));
273: PetscCall(PetscBLASIntCast(nmat,&nmat_));
274: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&nmat_,&k_,&sone,S,&lds_,fH+j*lda,&k_,&zero,h,&k_));
275: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,nmat,h,&M0));
276: PetscCall(BVSetActiveColumns(W,0,nmat));
277: PetscCall(BVMult(W,1.0,0.0,V,M0));
278: PetscCall(MatDestroy(&M0));
280: PetscCall(BVGetColumn(W,0,&w));
281: PetscCall(MatMult(A[0],w,Rv));
282: PetscCall(BVRestoreColumn(W,0,&w));
283: for (i=1;i<nmat;i++) {
284: PetscCall(BVGetColumn(W,i,&w));
285: PetscCall(MatMult(A[i],w,t));
286: PetscCall(BVRestoreColumn(W,i,&w));
287: PetscCall(VecAXPY(Rv,1.0,t));
288: }
289: /* Update right-hand side */
290: if (j) {
291: PetscCall(PetscBLASIntCast(ldh,&ldh_));
292: PetscCall(PetscArrayzero(Z,k*k));
293: PetscCall(PetscArrayzero(DS0,k*k));
294: PetscCall(PetscArraycpy(Z+(j-1)*k,dH+(j-1)*k,k));
295: /* Update DfH */
296: for (i=1;i<nmat;i++) {
297: if (i>1) {
298: beta = -g[i-1];
299: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,fH+(i-1)*k,&lda_,Z,&k_,&beta,DS0,&k_));
300: tt += -b[i-1];
301: for (ii=0;ii<k;ii++) H[ii+ii*ldh] += tt;
302: tt = b[i-1];
303: beta = 1.0/a[i-1];
304: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&beta,DS1,&k_,H,&ldh_,&beta,DS0,&k_));
305: F = DS0; DS0 = DS1; DS1 = F;
306: } else {
307: PetscCall(PetscArrayzero(DS1,k*k));
308: for (ii=0;ii<k;ii++) DS1[ii+(j-1)*k] = Z[ii+(j-1)*k]/a[0];
309: }
310: for (jj=j;jj<k;jj++) {
311: for (ii=0;ii<k;ii++) DfH[k*i+ii+jj*lda] += DS1[ii+jj*k];
312: }
313: }
314: for (ii=0;ii<k;ii++) H[ii+ii*ldh] += tt;
315: /* Update right-hand side */
316: PetscCall(PetscBLASIntCast(2*k,&k2_));
317: PetscCall(PetscBLASIntCast(j,&j_));
318: PetscCall(PetscBLASIntCast(k+rds,&krds_));
319: c0 = DS0;
320: PetscCall(PetscArrayzero(Rh,k));
321: for (i=0;i<nmat;i++) {
322: PetscCallBLAS("BLASgemv",BLASgemv_("N",&krds_,&j_,&sone,dVS,&k2_,fH+j*lda+i*k,&one,&zero,h,&one));
323: PetscCallBLAS("BLASgemv",BLASgemv_("N",&k_,&k_,&sone,S,&lds_,DfH+i*k+j*lda,&one,&sone,h,&one));
324: PetscCall(BVMultVec(V,1.0,0.0,t,h));
325: PetscCall(BVSetActiveColumns(dV,0,rds));
326: PetscCall(BVMultVec(dV,1.0,1.0,t,h+k));
327: PetscCall(BVGetColumn(W,i,&w));
328: PetscCall(MatMult(A[i],t,w));
329: PetscCall(BVRestoreColumn(W,i,&w));
330: if (i>0 && i<nmat-1) {
331: PetscCallBLAS("BLASgemv",BLASgemv_("C",&k_,&k_,&sone,S,&lds_,h,&one,&zero,c0,&one));
332: PetscCallBLAS("BLASgemv",BLASgemv_("C",&k_,&k_,&none,fH+i*k,&lda_,c0,&one,&sone,Rh,&one));
333: }
334: }
336: for (i=0;i<nmat;i++) h[i] = -1.0;
337: PetscCall(BVMultVec(W,1.0,1.0,Rv,h));
338: }
339: PetscCall(PetscFree4(h,DS0,DS1,Z));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static PetscErrorCode NRefSysSolve_mbe(PetscInt k,PetscInt sz,BV W,PetscScalar *w,BV Wt,PetscScalar *wt,PetscScalar *d,PetscScalar *dt,KSP ksp,BV T2,BV T3 ,PetscScalar *T4,PetscBool trans,Vec x1,PetscScalar *x2,Vec sol1,PetscScalar *sol2,Vec vw)
344: {
345: PetscInt i,j,incf,incc;
346: PetscScalar *y,*g,*xx2,*ww,y2,*dd;
347: Vec v,t,xx1;
348: BV WW,T;
350: PetscFunctionBegin;
351: PetscCall(PetscMalloc3(sz,&y,sz,&g,k,&xx2));
352: if (trans) {
353: WW = W; ww = w; dd = d; T = T3; incf = 0; incc = 1;
354: } else {
355: WW = Wt; ww = wt; dd = dt; T = T2; incf = 1; incc = 0;
356: }
357: xx1 = vw;
358: PetscCall(VecCopy(x1,xx1));
359: PetscCall(PetscArraycpy(xx2,x2,sz));
360: PetscCall(PetscArrayzero(sol2,k));
361: for (i=sz-1;i>=0;i--) {
362: PetscCall(BVGetColumn(WW,i,&v));
363: PetscCall(VecConjugate(v));
364: PetscCall(VecDot(xx1,v,y+i));
365: PetscCall(VecConjugate(v));
366: PetscCall(BVRestoreColumn(WW,i,&v));
367: for (j=0;j<i;j++) y[i] += ww[j+i*k]*xx2[j];
368: y[i] = -(y[i]-xx2[i])/dd[i];
369: PetscCall(BVGetColumn(T,i,&t));
370: PetscCall(VecAXPY(xx1,-y[i],t));
371: PetscCall(BVRestoreColumn(T,i,&t));
372: for (j=0;j<=i;j++) xx2[j] -= y[i]*T4[j*incf+incc*i+(i*incf+incc*j)*k];
373: g[i] = xx2[i];
374: }
375: if (trans) PetscCall(KSPSolveTranspose(ksp,xx1,sol1));
376: else PetscCall(KSPSolve(ksp,xx1,sol1));
377: if (trans) {
378: WW = Wt; ww = wt; dd = dt; T = T2; incf = 1; incc = 0;
379: } else {
380: WW = W; ww = w; dd = d; T = T3; incf = 0; incc = 1;
381: }
382: for (i=0;i<sz;i++) {
383: PetscCall(BVGetColumn(T,i,&t));
384: PetscCall(VecConjugate(t));
385: PetscCall(VecDot(sol1,t,&y2));
386: PetscCall(VecConjugate(t));
387: PetscCall(BVRestoreColumn(T,i,&t));
388: for (j=0;j<i;j++) y2 += sol2[j]*T4[j*incf+incc*i+(i*incf+incc*j)*k];
389: y2 = (g[i]-y2)/dd[i];
390: PetscCall(BVGetColumn(WW,i,&v));
391: PetscCall(VecAXPY(sol1,-y2,v));
392: for (j=0;j<i;j++) sol2[j] -= ww[j+i*k]*y2;
393: sol2[i] = y[i]+y2;
394: PetscCall(BVRestoreColumn(WW,i,&v));
395: }
396: PetscCall(PetscFree3(y,g,xx2));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode NRefSysSetup_mbe(PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,BV V,PEP_REFINE_EXPLICIT *matctx)
401: {
402: PetscInt i,j,l,nmat=pep->nmat,lda=nmat*k,deg=nmat-1;
403: Mat M1=matctx->M1,*A,*At,Mk;
404: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
405: PetscScalar s,ss,*DHii,*T12,*array,*Ts,*Tr,*M4=matctx->M4,sone=1.0,zero=0.0;
406: PetscScalar *w=matctx->w,*wt=matctx->wt,*d=matctx->d,*dt=matctx->dt;
407: PetscBLASInt lds_,lda_,k_;
408: MatStructure str;
409: PetscBool flg;
410: BV M2=matctx->M2,M3=matctx->M3,W=matctx->W,Wt=matctx->Wt;
411: Vec vc,vc2;
413: PetscFunctionBegin;
414: PetscCall(PetscMalloc3(nmat*k*k,&T12,k*k,&Tr,PetscMax(k*k,nmat),&Ts));
415: PetscCall(STGetMatStructure(pep->st,&str));
416: PetscCall(STGetTransform(pep->st,&flg));
417: if (flg) {
418: PetscCall(PetscMalloc1(pep->nmat,&At));
419: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
420: } else At = pep->A;
421: if (matctx->subc) A = matctx->A;
422: else A = At;
423: /* Form the explicit system matrix */
424: DHii = T12;
425: PetscCall(PetscArrayzero(DHii,k*k*nmat));
426: for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
427: for (l=2;l<nmat;l++) {
428: for (j=0;j<k;j++) {
429: for (i=0;i<k;i++) {
430: DHii[l*k+i+j*lda] = ((h-b[l-1])*DHii[(l-1)*k+i+j*lda]+fH[(l-1)*k+i+j*lda]-g[l-1]*DHii[(l-2)*k+i+j*lda])/a[l-1];
431: }
432: }
433: }
435: /* T11 */
436: if (!matctx->compM1) {
437: PetscCall(MatCopy(A[0],M1,DIFFERENT_NONZERO_PATTERN));
438: PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
439: for (j=1;j<nmat;j++) PetscCall(MatAXPY(M1,Ts[j],A[j],str));
440: }
442: /* T22 */
443: PetscCall(PetscBLASIntCast(lds,&lds_));
444: PetscCall(PetscBLASIntCast(k,&k_));
445: PetscCall(PetscBLASIntCast(lda,&lda_));
446: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
447: for (i=1;i<deg;i++) {
448: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
449: s = (i==1)?0.0:1.0;
450: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,M4,&k_));
451: }
453: /* T12 */
454: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&Mk));
455: for (i=1;i<nmat;i++) {
456: PetscCall(MatDenseGetArrayWrite(Mk,&array));
457: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,array,&k_));
458: PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
459: PetscCall(BVSetActiveColumns(W,0,k));
460: PetscCall(BVMult(W,1.0,0.0,V,Mk));
461: if (i==1) PetscCall(BVMatMult(W,A[i],M2));
462: else {
463: PetscCall(BVMatMult(W,A[i],M3)); /* using M3 as work space */
464: PetscCall(BVMult(M2,1.0,1.0,M3,NULL));
465: }
466: }
468: /* T21 */
469: PetscCall(MatDenseGetArrayWrite(Mk,&array));
470: for (i=1;i<deg;i++) {
471: s = (i==1)?0.0:1.0;
472: ss = PetscConj(fh[i]);
473: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,array,&k_));
474: }
475: PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
476: PetscCall(BVSetActiveColumns(M3,0,k));
477: PetscCall(BVMult(M3,1.0,0.0,V,Mk));
478: for (i=0;i<k;i++) {
479: PetscCall(BVGetColumn(M3,i,&vc));
480: PetscCall(VecConjugate(vc));
481: PetscCall(BVRestoreColumn(M3,i,&vc));
482: }
484: PetscCall(PEP_KSPSetOperators(ksp,M1,M1));
485: PetscCall(KSPSetUp(ksp));
486: PetscCall(MatDestroy(&Mk));
488: /* Set up for BEMW */
489: for (i=0;i<k;i++) {
490: PetscCall(BVGetColumn(M2,i,&vc));
491: PetscCall(BVGetColumn(W,i,&vc2));
492: PetscCall(NRefSysSolve_mbe(k,i,W,w,Wt,wt,d,dt,ksp,M2,M3,M4,PETSC_FALSE,vc,M4+i*k,vc2,w+i*k,matctx->t));
493: PetscCall(BVRestoreColumn(M2,i,&vc));
494: PetscCall(BVGetColumn(M3,i,&vc));
495: PetscCall(VecConjugate(vc));
496: PetscCall(VecDot(vc2,vc,&d[i]));
497: PetscCall(VecConjugate(vc));
498: PetscCall(BVRestoreColumn(M3,i,&vc));
499: for (j=0;j<i;j++) d[i] += M4[i+j*k]*w[j+i*k];
500: d[i] = M4[i+i*k]-d[i];
501: PetscCall(BVRestoreColumn(W,i,&vc2));
503: PetscCall(BVGetColumn(M3,i,&vc));
504: PetscCall(BVGetColumn(Wt,i,&vc2));
505: for (j=0;j<=i;j++) Ts[j] = M4[i+j*k];
506: PetscCall(NRefSysSolve_mbe(k,i,W,w,Wt,wt,d,dt,ksp,M2,M3,M4,PETSC_TRUE,vc,Ts,vc2,wt+i*k,matctx->t));
507: PetscCall(BVRestoreColumn(M3,i,&vc));
508: PetscCall(BVGetColumn(M2,i,&vc));
509: PetscCall(VecConjugate(vc2));
510: PetscCall(VecDot(vc,vc2,&dt[i]));
511: PetscCall(VecConjugate(vc2));
512: PetscCall(BVRestoreColumn(M2,i,&vc));
513: for (j=0;j<i;j++) dt[i] += M4[j+i*k]*wt[j+i*k];
514: dt[i] = M4[i+i*k]-dt[i];
515: PetscCall(BVRestoreColumn(Wt,i,&vc2));
516: }
518: if (flg) PetscCall(PetscFree(At));
519: PetscCall(PetscFree3(T12,Tr,Ts));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode NRefSysSetup_explicit(PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,BV V,PEP_REFINE_EXPLICIT *matctx,BV W)
524: {
525: PetscInt i,j,d,n,n0,m0,n1,m1,nmat=pep->nmat,lda=nmat*k,deg=nmat-1;
526: PetscInt *idxg=matctx->idxg,*idxp=matctx->idxp,idx,ncols;
527: Mat M,*E=matctx->E,*A,*At,Mk,Md;
528: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
529: PetscScalar s,ss,*DHii,*T22,*T21,*T12,*Ts,*Tr,*array,*ts,sone=1.0,zero=0.0;
530: PetscBLASInt lds_,lda_,k_;
531: const PetscInt *idxmc;
532: const PetscScalar *valsc,*carray;
533: MatStructure str;
534: Vec vc,vc0;
535: PetscBool flg;
537: PetscFunctionBegin;
538: PetscCall(PetscMalloc5(k*k,&T22,k*k,&T21,nmat*k*k,&T12,k*k,&Tr,k*k,&Ts));
539: PetscCall(STGetMatStructure(pep->st,&str));
540: PetscCall(KSPGetOperators(ksp,&M,NULL));
541: PetscCall(MatGetOwnershipRange(E[1],&n1,&m1));
542: PetscCall(MatGetOwnershipRange(E[0],&n0,&m0));
543: PetscCall(MatGetOwnershipRange(M,&n,NULL));
544: PetscCall(PetscMalloc1(nmat,&ts));
545: PetscCall(STGetTransform(pep->st,&flg));
546: if (flg) {
547: PetscCall(PetscMalloc1(pep->nmat,&At));
548: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
549: } else At = pep->A;
550: if (matctx->subc) A = matctx->A;
551: else A = At;
552: /* Form the explicit system matrix */
553: DHii = T12;
554: PetscCall(PetscArrayzero(DHii,k*k*nmat));
555: for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
556: for (d=2;d<nmat;d++) {
557: for (j=0;j<k;j++) {
558: for (i=0;i<k;i++) {
559: DHii[d*k+i+j*lda] = ((h-b[d-1])*DHii[(d-1)*k+i+j*lda]+fH[(d-1)*k+i+j*lda]-g[d-1]*DHii[(d-2)*k+i+j*lda])/a[d-1];
560: }
561: }
562: }
564: /* T11 */
565: if (!matctx->compM1) {
566: PetscCall(MatCopy(A[0],E[0],DIFFERENT_NONZERO_PATTERN));
567: PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
568: for (j=1;j<nmat;j++) PetscCall(MatAXPY(E[0],Ts[j],A[j],str));
569: }
570: for (i=n0;i<m0;i++) {
571: PetscCall(MatGetRow(E[0],i,&ncols,&idxmc,&valsc));
572: idx = n+i-n0;
573: for (j=0;j<ncols;j++) {
574: idxg[j] = matctx->map0[idxmc[j]];
575: }
576: PetscCall(MatSetValues(M,1,&idx,ncols,idxg,valsc,INSERT_VALUES));
577: PetscCall(MatRestoreRow(E[0],i,&ncols,&idxmc,&valsc));
578: }
580: /* T22 */
581: PetscCall(PetscBLASIntCast(lds,&lds_));
582: PetscCall(PetscBLASIntCast(k,&k_));
583: PetscCall(PetscBLASIntCast(lda,&lda_));
584: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
585: for (i=1;i<deg;i++) {
586: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
587: s = (i==1)?0.0:1.0;
588: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,T22,&k_));
589: }
590: for (j=0;j<k;j++) idxp[j] = matctx->map1[j];
591: for (i=0;i<m1-n1;i++) {
592: idx = n+m0-n0+i;
593: for (j=0;j<k;j++) {
594: Tr[j] = T22[n1+i+j*k];
595: }
596: PetscCall(MatSetValues(M,1,&idx,k,idxp,Tr,INSERT_VALUES));
597: }
599: /* T21 */
600: for (i=1;i<deg;i++) {
601: s = (i==1)?0.0:1.0;
602: ss = PetscConj(fh[i]);
603: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,T21,&k_));
604: }
605: PetscCall(BVSetActiveColumns(W,0,k));
606: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,T21,&Mk));
607: PetscCall(BVMult(W,1.0,0.0,V,Mk));
608: for (i=0;i<k;i++) {
609: PetscCall(BVGetColumn(W,i,&vc));
610: PetscCall(VecConjugate(vc));
611: PetscCall(VecGetArrayRead(vc,&carray));
612: idx = matctx->map1[i];
613: PetscCall(MatSetValues(M,1,&idx,m0-n0,matctx->map0+n0,carray,INSERT_VALUES));
614: PetscCall(VecRestoreArrayRead(vc,&carray));
615: PetscCall(BVRestoreColumn(W,i,&vc));
616: }
618: /* T12 */
619: for (i=1;i<nmat;i++) {
620: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,Ts,&k_));
621: for (j=0;j<k;j++) PetscCall(PetscArraycpy(T12+i*k+j*lda,Ts+j*k,k));
622: }
623: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,nmat-1,NULL,&Md));
624: for (i=0;i<nmat;i++) ts[i] = 1.0;
625: for (j=0;j<k;j++) {
626: PetscCall(MatDenseGetArrayWrite(Md,&array));
627: PetscCall(PetscArraycpy(array,T12+k+j*lda,(nmat-1)*k));
628: PetscCall(MatDenseRestoreArrayWrite(Md,&array));
629: PetscCall(BVSetActiveColumns(W,0,nmat-1));
630: PetscCall(BVMult(W,1.0,0.0,V,Md));
631: for (i=nmat-1;i>0;i--) {
632: PetscCall(BVGetColumn(W,i-1,&vc0));
633: PetscCall(BVGetColumn(W,i,&vc));
634: PetscCall(MatMult(A[i],vc0,vc));
635: PetscCall(BVRestoreColumn(W,i-1,&vc0));
636: PetscCall(BVRestoreColumn(W,i,&vc));
637: }
638: PetscCall(BVSetActiveColumns(W,1,nmat));
639: PetscCall(BVGetColumn(W,0,&vc0));
640: PetscCall(BVMultVec(W,1.0,0.0,vc0,ts));
641: PetscCall(VecGetArrayRead(vc0,&carray));
642: idx = matctx->map1[j];
643: PetscCall(MatSetValues(M,m0-n0,matctx->map0+n0,1,&idx,carray,INSERT_VALUES));
644: PetscCall(VecRestoreArrayRead(vc0,&carray));
645: PetscCall(BVRestoreColumn(W,0,&vc0));
646: }
647: PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
648: PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
649: PetscCall(PEP_KSPSetOperators(ksp,M,M));
650: PetscCall(KSPSetUp(ksp));
651: PetscCall(PetscFree(ts));
652: PetscCall(MatDestroy(&Mk));
653: PetscCall(MatDestroy(&Md));
654: if (flg) PetscCall(PetscFree(At));
655: PetscCall(PetscFree5(T22,T21,T12,Tr,Ts));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode NRefSysSolve_explicit(PetscInt k,KSP ksp,Vec Rv,PetscScalar *Rh,Vec dVi,PetscScalar *dHi,PEP_REFINE_EXPLICIT *matctx)
660: {
661: PetscInt n0,m0,n1,m1,i;
662: PetscScalar *arrayV;
663: const PetscScalar *array;
665: PetscFunctionBegin;
666: PetscCall(MatGetOwnershipRange(matctx->E[1],&n1,&m1));
667: PetscCall(MatGetOwnershipRange(matctx->E[0],&n0,&m0));
669: /* Right side */
670: PetscCall(VecGetArrayRead(Rv,&array));
671: PetscCall(VecSetValues(matctx->tN,m0-n0,matctx->map0+n0,array,INSERT_VALUES));
672: PetscCall(VecRestoreArrayRead(Rv,&array));
673: PetscCall(VecSetValues(matctx->tN,m1-n1,matctx->map1+n1,Rh+n1,INSERT_VALUES));
674: PetscCall(VecAssemblyBegin(matctx->tN));
675: PetscCall(VecAssemblyEnd(matctx->tN));
677: /* Solve */
678: PetscCall(KSPSolve(ksp,matctx->tN,matctx->ttN));
680: /* Retrieve solution */
681: PetscCall(VecGetArray(dVi,&arrayV));
682: PetscCall(VecGetArrayRead(matctx->ttN,&array));
683: PetscCall(PetscArraycpy(arrayV,array,m0-n0));
684: PetscCall(VecRestoreArray(dVi,&arrayV));
685: if (!matctx->subc) {
686: PetscCall(VecGetArray(matctx->t1,&arrayV));
687: for (i=0;i<m1-n1;i++) arrayV[i] = array[m0-n0+i];
688: PetscCall(VecRestoreArray(matctx->t1,&arrayV));
689: PetscCall(VecRestoreArrayRead(matctx->ttN,&array));
690: PetscCall(VecScatterBegin(matctx->scatterctx,matctx->t1,matctx->vseq,INSERT_VALUES,SCATTER_FORWARD));
691: PetscCall(VecScatterEnd(matctx->scatterctx,matctx->t1,matctx->vseq,INSERT_VALUES,SCATTER_FORWARD));
692: PetscCall(VecGetArrayRead(matctx->vseq,&array));
693: for (i=0;i<k;i++) dHi[i] = array[i];
694: PetscCall(VecRestoreArrayRead(matctx->vseq,&array));
695: }
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: static PetscErrorCode NRefSysIter(PetscInt i,PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar *H,PetscInt ldh,Vec Rv,PetscScalar *Rh,BV V,Vec dVi,PetscScalar *dHi,PEP_REFINE_EXPLICIT *matctx,BV W)
700: {
701: PetscInt j,m,lda=pep->nmat*k,n0,m0,idx;
702: PetscMPIInt root,len;
703: PetscScalar *array2,h;
704: const PetscScalar *array;
705: Vec R,Vi;
706: PEP_REFINE_MATSHELL *ctx;
707: Mat M;
709: PetscFunctionBegin;
710: if (!matctx || !matctx->subc) {
711: for (j=0;j<pep->nmat;j++) fh[j] = fH[j*k+i+i*lda];
712: h = H[i+i*ldh];
713: idx = i;
714: R = Rv;
715: Vi = dVi;
716: switch (pep->scheme) {
717: case PEP_REFINE_SCHEME_EXPLICIT:
718: PetscCall(NRefSysSetup_explicit(pep,k,ksp,fH,S,lds,fh,h,V,matctx,W));
719: matctx->compM1 = PETSC_FALSE;
720: break;
721: case PEP_REFINE_SCHEME_MBE:
722: PetscCall(NRefSysSetup_mbe(pep,k,ksp,fH,S,lds,fh,h,V,matctx));
723: matctx->compM1 = PETSC_FALSE;
724: break;
725: case PEP_REFINE_SCHEME_SCHUR:
726: PetscCall(KSPGetOperators(ksp,&M,NULL));
727: PetscCall(MatShellGetContext(M,&ctx));
728: PetscCall(NRefSysSetup_shell(pep,k,fH,S,lds,fh,h,ctx));
729: ctx->compM1 = PETSC_FALSE;
730: break;
731: }
732: } else {
733: if (i%matctx->subc->n==0 && (idx=i+matctx->subc->color)<k) {
734: for (j=0;j<pep->nmat;j++) fh[j] = fH[j*k+idx+idx*lda];
735: h = H[idx+idx*ldh];
736: matctx->idx = idx;
737: switch (pep->scheme) {
738: case PEP_REFINE_SCHEME_EXPLICIT:
739: PetscCall(NRefSysSetup_explicit(pep,k,ksp,fH,S,lds,fh,h,matctx->V,matctx,matctx->W));
740: matctx->compM1 = PETSC_FALSE;
741: break;
742: case PEP_REFINE_SCHEME_MBE:
743: PetscCall(NRefSysSetup_mbe(pep,k,ksp,fH,S,lds,fh,h,matctx->V,matctx));
744: matctx->compM1 = PETSC_FALSE;
745: break;
746: case PEP_REFINE_SCHEME_SCHUR:
747: break;
748: }
749: } else idx = matctx->idx;
750: PetscCall(VecScatterBegin(matctx->scatter_id[i%matctx->subc->n],Rv,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
751: PetscCall(VecScatterEnd(matctx->scatter_id[i%matctx->subc->n],Rv,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
752: PetscCall(VecGetArrayRead(matctx->tg,&array));
753: PetscCall(VecPlaceArray(matctx->t,array));
754: PetscCall(VecCopy(matctx->t,matctx->Rv));
755: PetscCall(VecResetArray(matctx->t));
756: PetscCall(VecRestoreArrayRead(matctx->tg,&array));
757: R = matctx->Rv;
758: Vi = matctx->Vi;
759: }
760: if (idx==i && idx<k) {
761: switch (pep->scheme) {
762: case PEP_REFINE_SCHEME_EXPLICIT:
763: PetscCall(NRefSysSolve_explicit(k,ksp,R,Rh,Vi,dHi,matctx));
764: break;
765: case PEP_REFINE_SCHEME_MBE:
766: PetscCall(NRefSysSolve_mbe(k,k,matctx->W,matctx->w,matctx->Wt,matctx->wt,matctx->d,matctx->dt,ksp,matctx->M2,matctx->M3 ,matctx->M4,PETSC_FALSE,R,Rh,Vi,dHi,matctx->t));
767: break;
768: case PEP_REFINE_SCHEME_SCHUR:
769: PetscCall(NRefSysSolve_shell(ksp,pep->nmat,R,Rh,k,Vi,dHi));
770: break;
771: }
772: }
773: if (matctx && matctx->subc) {
774: PetscCall(VecGetLocalSize(Vi,&m));
775: PetscCall(VecGetArrayRead(Vi,&array));
776: PetscCall(VecGetArray(matctx->tg,&array2));
777: PetscCall(PetscArraycpy(array2,array,m));
778: PetscCall(VecRestoreArray(matctx->tg,&array2));
779: PetscCall(VecRestoreArrayRead(Vi,&array));
780: PetscCall(VecScatterBegin(matctx->scatter_id[i%matctx->subc->n],matctx->tg,dVi,INSERT_VALUES,SCATTER_REVERSE));
781: PetscCall(VecScatterEnd(matctx->scatter_id[i%matctx->subc->n],matctx->tg,dVi,INSERT_VALUES,SCATTER_REVERSE));
782: switch (pep->scheme) {
783: case PEP_REFINE_SCHEME_EXPLICIT:
784: PetscCall(MatGetOwnershipRange(matctx->E[0],&n0,&m0));
785: PetscCall(VecGetArrayRead(matctx->ttN,&array));
786: PetscCall(VecPlaceArray(matctx->tp,array+m0-n0));
787: PetscCall(VecScatterBegin(matctx->scatterp_id[i%matctx->subc->n],matctx->tp,matctx->tpg,INSERT_VALUES,SCATTER_FORWARD));
788: PetscCall(VecScatterEnd(matctx->scatterp_id[i%matctx->subc->n],matctx->tp,matctx->tpg,INSERT_VALUES,SCATTER_FORWARD));
789: PetscCall(VecResetArray(matctx->tp));
790: PetscCall(VecRestoreArrayRead(matctx->ttN,&array));
791: PetscCall(VecGetArrayRead(matctx->tpg,&array));
792: for (j=0;j<k;j++) dHi[j] = array[j];
793: PetscCall(VecRestoreArrayRead(matctx->tpg,&array));
794: break;
795: case PEP_REFINE_SCHEME_MBE:
796: root = 0;
797: for (j=0;j<i%matctx->subc->n;j++) root += matctx->subc->subsize[j];
798: PetscCall(PetscMPIIntCast(k,&len));
799: PetscCallMPI(MPI_Bcast(dHi,len,MPIU_SCALAR,root,matctx->subc->dupparent));
800: break;
801: case PEP_REFINE_SCHEME_SCHUR:
802: break;
803: }
804: }
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: static PetscErrorCode PEPNRefForwardSubstitution(PEP pep,PetscInt k,PetscScalar *S,PetscInt lds,PetscScalar *H,PetscInt ldh,PetscScalar *fH,BV dV,PetscScalar *dVS,PetscInt *rds,PetscScalar *dH,PetscInt lddh,KSP ksp,PEP_REFINE_EXPLICIT *matctx)
809: {
810: PetscInt i,nmat=pep->nmat,lda=nmat*k;
811: PetscScalar *fh,*Rh,*DfH;
812: PetscReal norm;
813: BV W;
814: Vec Rv,t,dvi;
815: PEP_REFINE_MATSHELL *ctx;
816: Mat M,*At;
817: PetscBool flg,lindep;
819: PetscFunctionBegin;
820: PetscCall(PetscMalloc2(nmat*k*k,&DfH,k,&Rh));
821: *rds = 0;
822: PetscCall(BVCreateVec(pep->V,&Rv));
823: switch (pep->scheme) {
824: case PEP_REFINE_SCHEME_EXPLICIT:
825: PetscCall(BVCreateVec(pep->V,&t));
826: PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
827: PetscCall(PetscMalloc1(nmat,&fh));
828: break;
829: case PEP_REFINE_SCHEME_MBE:
830: if (matctx->subc) {
831: PetscCall(BVCreateVec(pep->V,&t));
832: PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
833: } else {
834: W = matctx->W;
835: PetscCall(PetscObjectReference((PetscObject)W));
836: t = matctx->t;
837: PetscCall(PetscObjectReference((PetscObject)t));
838: }
839: PetscCall(BVScale(matctx->W,0.0));
840: PetscCall(BVScale(matctx->Wt,0.0));
841: PetscCall(BVScale(matctx->M2,0.0));
842: PetscCall(BVScale(matctx->M3,0.0));
843: PetscCall(PetscMalloc1(nmat,&fh));
844: break;
845: case PEP_REFINE_SCHEME_SCHUR:
846: PetscCall(KSPGetOperators(ksp,&M,NULL));
847: PetscCall(MatShellGetContext(M,&ctx));
848: PetscCall(BVCreateVec(pep->V,&t));
849: PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
850: fh = ctx->fih;
851: break;
852: }
853: PetscCall(PetscArrayzero(dVS,2*k*k));
854: PetscCall(PetscArrayzero(DfH,lda*k));
855: PetscCall(STGetTransform(pep->st,&flg));
856: if (flg) {
857: PetscCall(PetscMalloc1(pep->nmat,&At));
858: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
859: } else At = pep->A;
861: /* Main loop for computing the i-th columns of dX and dS */
862: for (i=0;i<k;i++) {
863: /* Compute and update i-th column of the right hand side */
864: PetscCall(PetscArrayzero(Rh,k));
865: PetscCall(NRefRightSide(nmat,pep->pbc,At,k,pep->V,S,lds,i,H,ldh,fH,DfH,dH,dV,dVS,*rds,Rv,Rh,W,t));
867: /* Update and solve system */
868: PetscCall(BVGetColumn(dV,i,&dvi));
869: PetscCall(NRefSysIter(i,pep,k,ksp,fH,S,lds,fh,H,ldh,Rv,Rh,pep->V,dvi,dH+i*k,matctx,W));
870: /* Orthogonalize computed solution */
871: PetscCall(BVOrthogonalizeVec(pep->V,dvi,dVS+i*2*k,&norm,&lindep));
872: PetscCall(BVRestoreColumn(dV,i,&dvi));
873: if (!lindep) {
874: PetscCall(BVOrthogonalizeColumn(dV,i,dVS+k+i*2*k,&norm,&lindep));
875: if (!lindep) {
876: dVS[k+i+i*2*k] = norm;
877: PetscCall(BVScaleColumn(dV,i,1.0/norm));
878: (*rds)++;
879: }
880: }
881: }
882: PetscCall(BVSetActiveColumns(dV,0,*rds));
883: PetscCall(VecDestroy(&t));
884: PetscCall(VecDestroy(&Rv));
885: PetscCall(BVDestroy(&W));
886: if (flg) PetscCall(PetscFree(At));
887: PetscCall(PetscFree2(DfH,Rh));
888: if (pep->scheme!=PEP_REFINE_SCHEME_SCHUR) PetscCall(PetscFree(fh));
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: static PetscErrorCode NRefOrthogStep(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *S,PetscInt lds)
893: {
894: PetscInt j,nmat=pep->nmat,deg=nmat-1,lda=nmat*k,ldg;
895: PetscScalar *G,*tau,sone=1.0,zero=0.0,*work;
896: PetscBLASInt lds_,k_,ldh_,ldg_,lda_;
898: PetscFunctionBegin;
899: PetscCall(PetscMalloc3(k,&tau,k,&work,deg*k*k,&G));
900: PetscCall(PetscBLASIntCast(lds,&lds_));
901: PetscCall(PetscBLASIntCast(lda,&lda_));
902: PetscCall(PetscBLASIntCast(k,&k_));
904: /* Form auxiliary matrix for the orthogonalization step */
905: ldg = deg*k;
906: PetscCall(PEPEvaluateBasisforMatrix(pep,nmat,k,H,ldh,fH));
907: PetscCall(PetscBLASIntCast(ldg,&ldg_));
908: PetscCall(PetscBLASIntCast(ldh,&ldh_));
909: for (j=0;j<deg;j++) {
910: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,fH+j*k,&lda_,&zero,G+j*k,&ldg_));
911: }
912: /* Orthogonalize and update S */
913: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
914: PetscCallLAPACKInfo("LAPACKgeqrf",LAPACKgeqrf_(&ldg_,&k_,G,&ldg_,tau,work,&k_,&info));
915: PetscCall(PetscFPTrapPop());
916: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&k_,&k_,&sone,G,&ldg_,S,&lds_));
918: /* Update H */
919: PetscCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&k_,&k_,&sone,G,&ldg_,H,&ldh_));
920: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&k_,&k_,&sone,G,&ldg_,H,&ldh_));
921: PetscCall(PetscFree3(tau,work,G));
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: static PetscErrorCode PEPNRefUpdateInvPair(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *dH,PetscScalar *S,PetscInt lds,BV dV,PetscScalar *dVS,PetscInt rds)
926: {
927: PetscInt i,j,nmat=pep->nmat,lda=nmat*k;
928: PetscScalar *tau,*array,*work;
929: PetscBLASInt lds_,k_,lda_,ldh_,kdrs_,k2_;
930: Mat M0;
932: PetscFunctionBegin;
933: PetscCall(PetscMalloc2(k,&tau,k,&work));
934: PetscCall(PetscBLASIntCast(lds,&lds_));
935: PetscCall(PetscBLASIntCast(lda,&lda_));
936: PetscCall(PetscBLASIntCast(ldh,&ldh_));
937: PetscCall(PetscBLASIntCast(k,&k_));
938: PetscCall(PetscBLASIntCast(2*k,&k2_));
939: PetscCall(PetscBLASIntCast((k+rds),&kdrs_));
940: /* Update H */
941: for (j=0;j<k;j++) {
942: for (i=0;i<k;i++) H[i+j*ldh] -= dH[i+j*k];
943: }
944: /* Update V */
945: for (j=0;j<k;j++) {
946: for (i=0;i<k;i++) dVS[i+j*2*k] = -dVS[i+j*2*k]+S[i+j*lds];
947: for (i=k;i<2*k;i++) dVS[i+j*2*k] = -dVS[i+j*2*k];
948: }
949: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
950: PetscCallLAPACKInfo("LAPACKgeqrf",LAPACKgeqrf_(&kdrs_,&k_,dVS,&k2_,tau,work,&k_,&info));
951: /* Copy triangular matrix in S */
952: for (j=0;j<k;j++) {
953: for (i=0;i<=j;i++) S[i+j*lds] = dVS[i+j*2*k];
954: for (i=j+1;i<k;i++) S[i+j*lds] = 0.0;
955: }
956: PetscCallLAPACKInfo("LAPACKorgqr",LAPACKorgqr_(&k2_,&k_,&k_,dVS,&k2_,tau,work,&k_,&info));
957: PetscCall(PetscFPTrapPop());
958: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M0));
959: PetscCall(MatDenseGetArrayWrite(M0,&array));
960: for (j=0;j<k;j++) PetscCall(PetscArraycpy(array+j*k,dVS+j*2*k,k));
961: PetscCall(MatDenseRestoreArrayWrite(M0,&array));
962: PetscCall(BVMultInPlace(pep->V,M0,0,k));
963: if (rds) {
964: PetscCall(MatDenseGetArrayWrite(M0,&array));
965: for (j=0;j<k;j++) PetscCall(PetscArraycpy(array+j*k,dVS+k+j*2*k,rds));
966: PetscCall(MatDenseRestoreArrayWrite(M0,&array));
967: PetscCall(BVMultInPlace(dV,M0,0,k));
968: PetscCall(BVMult(pep->V,1.0,1.0,dV,NULL));
969: }
970: PetscCall(MatDestroy(&M0));
971: PetscCall(NRefOrthogStep(pep,k,H,ldh,fH,S,lds));
972: PetscCall(PetscFree2(tau,work));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: static PetscErrorCode PEPNRefSetUp(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PEP_REFINE_EXPLICIT *matctx,PetscBool ini)
977: {
978: PEP_REFINE_MATSHELL *ctx;
979: PetscScalar t,*coef;
980: const PetscScalar *array;
981: MatStructure str;
982: PetscInt j,nmat=pep->nmat,n0,m0,n1,m1,n0_,m0_,n1_,m1_,N0,N1,p,*idx1,*idx2,count,si,i,l0;
983: MPI_Comm comm;
984: PetscMPIInt np;
985: const PetscInt *rgs0,*rgs1;
986: Mat B,C,*E,*A,*At;
987: IS is1,is2;
988: Vec v;
989: PetscBool flg;
990: Mat M,P;
992: PetscFunctionBegin;
993: PetscCall(PetscMalloc1(nmat,&coef));
994: PetscCall(STGetTransform(pep->st,&flg));
995: if (flg) {
996: PetscCall(PetscMalloc1(pep->nmat,&At));
997: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
998: } else At = pep->A;
999: switch (pep->scheme) {
1000: case PEP_REFINE_SCHEME_EXPLICIT:
1001: if (ini) {
1002: if (matctx->subc) {
1003: A = matctx->A;
1004: PetscCall(PetscSubcommGetChild(matctx->subc,&comm));
1005: } else {
1006: A = At;
1007: PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
1008: }
1009: E = matctx->E;
1010: PetscCall(STGetMatStructure(pep->st,&str));
1011: PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&E[0]));
1012: j = matctx->subc?matctx->subc->color:0;
1013: PetscCall(PEPEvaluateBasis(pep,H[j+j*ldh],0,coef,NULL));
1014: for (j=1;j<nmat;j++) PetscCall(MatAXPY(E[0],coef[j],A[j],str));
1015: PetscCall(MatCreateDense(comm,PETSC_DECIDE,PETSC_DECIDE,k,k,NULL,&E[1]));
1016: PetscCall(MatGetOwnershipRange(E[0],&n0,&m0));
1017: PetscCall(MatGetOwnershipRange(E[1],&n1,&m1));
1018: PetscCall(MatGetOwnershipRangeColumn(E[0],&n0_,&m0_));
1019: PetscCall(MatGetOwnershipRangeColumn(E[1],&n1_,&m1_));
1020: /* T12 and T21 are computed from V and V*, so,
1021: they must have the same column and row ranges */
1022: PetscCheck(m0_-n0_==m0-n0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent dimensions");
1023: PetscCall(MatCreateDense(comm,m0-n0,m1_-n1_,PETSC_DECIDE,PETSC_DECIDE,NULL,&B));
1024: PetscCall(MatCreateDense(comm,m1-n1,m0_-n0_,PETSC_DECIDE,PETSC_DECIDE,NULL,&C));
1025: PetscCall(MatCreateTile(1.0,E[0],1.0,B,1.0,C,1.0,E[1],&M));
1026: PetscCall(MatDestroy(&B));
1027: PetscCall(MatDestroy(&C));
1028: matctx->compM1 = PETSC_TRUE;
1029: PetscCall(MatGetSize(E[0],NULL,&N0));
1030: PetscCall(MatGetSize(E[1],NULL,&N1));
1031: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)M),&np));
1032: PetscCall(MatGetOwnershipRanges(E[0],&rgs0));
1033: PetscCall(MatGetOwnershipRanges(E[1],&rgs1));
1034: PetscCall(PetscMalloc4(PetscMax(k,N1),&matctx->idxp,N0,&matctx->idxg,N0,&matctx->map0,N1,&matctx->map1));
1035: /* Create column (and row) mapping */
1036: for (p=0;p<np;p++) {
1037: for (j=rgs0[p];j<rgs0[p+1];j++) matctx->map0[j] = j+rgs1[p];
1038: for (j=rgs1[p];j<rgs1[p+1];j++) matctx->map1[j] = j+rgs0[p+1];
1039: }
1040: PetscCall(MatCreateVecs(M,NULL,&matctx->tN));
1041: PetscCall(MatCreateVecs(matctx->E[1],NULL,&matctx->t1));
1042: PetscCall(VecDuplicate(matctx->tN,&matctx->ttN));
1043: if (matctx->subc) {
1044: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
1045: count = np*k;
1046: PetscCall(PetscMalloc2(count,&idx1,count,&idx2));
1047: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),m1-n1,PETSC_DECIDE,&matctx->tp));
1048: PetscCall(VecGetOwnershipRange(matctx->tp,&l0,NULL));
1049: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),k,PETSC_DECIDE,&matctx->tpg));
1050: for (si=0;si<matctx->subc->n;si++) {
1051: if (matctx->subc->color==si) {
1052: j=0;
1053: if (matctx->subc->color==si) {
1054: for (p=0;p<np;p++) {
1055: for (i=n1;i<m1;i++) {
1056: idx1[j] = l0+i-n1;
1057: idx2[j++] =p*k+i;
1058: }
1059: }
1060: }
1061: count = np*(m1-n1);
1062: } else count =0;
1063: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),count,idx1,PETSC_COPY_VALUES,&is1));
1064: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),count,idx2,PETSC_COPY_VALUES,&is2));
1065: PetscCall(VecScatterCreate(matctx->tp,is1,matctx->tpg,is2,&matctx->scatterp_id[si]));
1066: PetscCall(ISDestroy(&is1));
1067: PetscCall(ISDestroy(&is2));
1068: }
1069: PetscCall(PetscFree2(idx1,idx2));
1070: } else PetscCall(VecScatterCreateToAll(matctx->t1,&matctx->scatterctx,&matctx->vseq));
1071: P = M;
1072: } else {
1073: if (matctx->subc) {
1074: /* Scatter vectors pep->V */
1075: for (i=0;i<k;i++) {
1076: PetscCall(BVGetColumn(pep->V,i,&v));
1077: PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1078: PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1079: PetscCall(BVRestoreColumn(pep->V,i,&v));
1080: PetscCall(VecGetArrayRead(matctx->tg,&array));
1081: PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
1082: PetscCall(BVInsertVec(matctx->V,i,matctx->t));
1083: PetscCall(VecResetArray(matctx->t));
1084: PetscCall(VecRestoreArrayRead(matctx->tg,&array));
1085: }
1086: }
1087: }
1088: break;
1089: case PEP_REFINE_SCHEME_MBE:
1090: if (ini) {
1091: if (matctx->subc) {
1092: A = matctx->A;
1093: PetscCall(PetscSubcommGetChild(matctx->subc,&comm));
1094: } else {
1095: matctx->V = pep->V;
1096: A = At;
1097: PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
1098: PetscCall(MatCreateVecs(pep->A[0],&matctx->t,NULL));
1099: }
1100: PetscCall(STGetMatStructure(pep->st,&str));
1101: PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&matctx->M1));
1102: j = matctx->subc?matctx->subc->color:0;
1103: PetscCall(PEPEvaluateBasis(pep,H[j+j*ldh],0,coef,NULL));
1104: for (j=1;j<nmat;j++) PetscCall(MatAXPY(matctx->M1,coef[j],A[j],str));
1105: PetscCall(BVDuplicateResize(matctx->V,PetscMax(k,pep->nmat),&matctx->W));
1106: PetscCall(BVDuplicateResize(matctx->V,k,&matctx->M2));
1107: PetscCall(BVDuplicate(matctx->M2,&matctx->M3));
1108: PetscCall(BVDuplicate(matctx->M2,&matctx->Wt));
1109: PetscCall(PetscMalloc5(k*k,&matctx->M4,k*k,&matctx->w,k*k,&matctx->wt,k,&matctx->d,k,&matctx->dt));
1110: matctx->compM1 = PETSC_TRUE;
1111: M = matctx->M1;
1112: P = M;
1113: }
1114: break;
1115: case PEP_REFINE_SCHEME_SCHUR:
1116: if (ini) {
1117: PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
1118: PetscCall(MatGetSize(At[0],&m0,&n0));
1119: PetscCall(PetscMalloc1(1,&ctx));
1120: PetscCall(STGetMatStructure(pep->st,&str));
1121: /* Create a shell matrix to solve the linear system */
1122: ctx->V = pep->V;
1123: ctx->k = k; ctx->nmat = nmat;
1124: PetscCall(PetscMalloc5(nmat,&ctx->A,k*k,&ctx->M4,k,&ctx->pM4,2*k*k,&ctx->work,nmat,&ctx->fih));
1125: for (i=0;i<nmat;i++) ctx->A[i] = At[i];
1126: PetscCall(PetscArrayzero(ctx->M4,k*k));
1127: PetscCall(MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,m0,n0,ctx,&M));
1128: PetscCall(MatShellSetOperation(M,MATOP_MULT,(PetscErrorCodeFn*)MatMult_FS));
1129: PetscCall(BVDuplicateResize(ctx->V,PetscMax(k,pep->nmat),&ctx->W));
1130: PetscCall(BVDuplicateResize(ctx->V,k,&ctx->M2));
1131: PetscCall(BVDuplicate(ctx->M2,&ctx->M3));
1132: PetscCall(BVCreateVec(pep->V,&ctx->t));
1133: PetscCall(MatDuplicate(At[0],MAT_COPY_VALUES,&ctx->M1));
1134: PetscCall(PEPEvaluateBasis(pep,H[0],0,coef,NULL));
1135: for (j=1;j<nmat;j++) PetscCall(MatAXPY(ctx->M1,coef[j],At[j],str));
1136: PetscCall(MatDuplicate(At[0],MAT_COPY_VALUES,&P));
1137: /* Compute a precond matrix for the system */
1138: t = H[0];
1139: PetscCall(PEPEvaluateBasis(pep,t,0,coef,NULL));
1140: for (j=1;j<nmat;j++) PetscCall(MatAXPY(P,coef[j],At[j],str));
1141: ctx->compM1 = PETSC_TRUE;
1142: }
1143: break;
1144: }
1145: if (ini) {
1146: PetscCall(PEPRefineGetKSP(pep,&pep->refineksp));
1147: PetscCall(KSPSetErrorIfNotConverged(pep->refineksp,PETSC_TRUE));
1148: PetscCall(PEP_KSPSetOperators(pep->refineksp,M,P));
1149: PetscCall(KSPSetFromOptions(pep->refineksp));
1150: }
1152: if (!ini && matctx && matctx->subc) {
1153: /* Scatter vectors pep->V */
1154: for (i=0;i<k;i++) {
1155: PetscCall(BVGetColumn(pep->V,i,&v));
1156: PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1157: PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1158: PetscCall(BVRestoreColumn(pep->V,i,&v));
1159: PetscCall(VecGetArrayRead(matctx->tg,&array));
1160: PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
1161: PetscCall(BVInsertVec(matctx->V,i,matctx->t));
1162: PetscCall(VecResetArray(matctx->t));
1163: PetscCall(VecRestoreArrayRead(matctx->tg,&array));
1164: }
1165: }
1166: PetscCall(PetscFree(coef));
1167: if (flg) PetscCall(PetscFree(At));
1168: PetscFunctionReturn(PETSC_SUCCESS);
1169: }
1171: static PetscErrorCode NRefSubcommSetup(PEP pep,PetscInt k,PEP_REFINE_EXPLICIT *matctx,PetscInt nsubc)
1172: {
1173: PetscInt i,si,j,m0,n0,nloc0,nloc_sub,*idx1,*idx2;
1174: IS is1,is2;
1175: BVType type;
1176: Vec v;
1177: const PetscScalar *array;
1178: Mat *A;
1179: PetscBool flg;
1180: MPI_Comm contpar,child;
1182: PetscFunctionBegin;
1183: PetscCall(STGetTransform(pep->st,&flg));
1184: if (flg) {
1185: PetscCall(PetscMalloc1(pep->nmat,&A));
1186: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&A[i]));
1187: } else A = pep->A;
1188: PetscCall(PetscSubcommGetChild(matctx->subc,&child));
1189: PetscCall(PetscSubcommGetContiguousParent(matctx->subc,&contpar));
1191: /* Duplicate pep matrices */
1192: PetscCall(PetscMalloc3(pep->nmat,&matctx->A,nsubc,&matctx->scatter_id,nsubc,&matctx->scatterp_id));
1193: for (i=0;i<pep->nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],0,child,MAT_INITIAL_MATRIX,&matctx->A[i]));
1195: /* Create Scatter */
1196: PetscCall(MatCreateVecs(matctx->A[0],&matctx->t,NULL));
1197: PetscCall(MatGetLocalSize(matctx->A[0],&nloc_sub,NULL));
1198: PetscCall(VecCreateMPI(contpar,nloc_sub,PETSC_DECIDE,&matctx->tg));
1199: PetscCall(BVGetColumn(pep->V,0,&v));
1200: PetscCall(VecGetOwnershipRange(v,&n0,&m0));
1201: nloc0 = m0-n0;
1202: PetscCall(PetscMalloc2(matctx->subc->n*nloc0,&idx1,matctx->subc->n*nloc0,&idx2));
1203: j = 0;
1204: for (si=0;si<matctx->subc->n;si++) {
1205: for (i=n0;i<m0;i++) {
1206: idx1[j] = i;
1207: idx2[j++] = i+pep->n*si;
1208: }
1209: }
1210: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),matctx->subc->n*nloc0,idx1,PETSC_COPY_VALUES,&is1));
1211: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),matctx->subc->n*nloc0,idx2,PETSC_COPY_VALUES,&is2));
1212: PetscCall(VecScatterCreate(v,is1,matctx->tg,is2,&matctx->scatter_sub));
1213: PetscCall(ISDestroy(&is1));
1214: PetscCall(ISDestroy(&is2));
1215: for (si=0;si<matctx->subc->n;si++) {
1216: j=0;
1217: for (i=n0;i<m0;i++) {
1218: idx1[j] = i;
1219: idx2[j++] = i+pep->n*si;
1220: }
1221: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),nloc0,idx1,PETSC_COPY_VALUES,&is1));
1222: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),nloc0,idx2,PETSC_COPY_VALUES,&is2));
1223: PetscCall(VecScatterCreate(v,is1,matctx->tg,is2,&matctx->scatter_id[si]));
1224: PetscCall(ISDestroy(&is1));
1225: PetscCall(ISDestroy(&is2));
1226: }
1227: PetscCall(BVRestoreColumn(pep->V,0,&v));
1228: PetscCall(PetscFree2(idx1,idx2));
1230: /* Duplicate pep->V vecs */
1231: PetscCall(BVGetType(pep->V,&type));
1232: PetscCall(BVCreate(child,&matctx->V));
1233: PetscCall(BVSetType(matctx->V,type));
1234: PetscCall(BVSetSizesFromVec(matctx->V,matctx->t,k));
1235: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) PetscCall(BVDuplicateResize(matctx->V,PetscMax(k,pep->nmat),&matctx->W));
1236: for (i=0;i<k;i++) {
1237: PetscCall(BVGetColumn(pep->V,i,&v));
1238: PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1239: PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
1240: PetscCall(BVRestoreColumn(pep->V,i,&v));
1241: PetscCall(VecGetArrayRead(matctx->tg,&array));
1242: PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
1243: PetscCall(BVInsertVec(matctx->V,i,matctx->t));
1244: PetscCall(VecResetArray(matctx->t));
1245: PetscCall(VecRestoreArrayRead(matctx->tg,&array));
1246: }
1248: PetscCall(VecDuplicate(matctx->t,&matctx->Rv));
1249: PetscCall(VecDuplicate(matctx->t,&matctx->Vi));
1250: if (flg) PetscCall(PetscFree(A));
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: static PetscErrorCode NRefSubcommDestroy(PEP pep,PEP_REFINE_EXPLICIT *matctx)
1255: {
1256: PetscInt i;
1258: PetscFunctionBegin;
1259: PetscCall(VecScatterDestroy(&matctx->scatter_sub));
1260: for (i=0;i<matctx->subc->n;i++) PetscCall(VecScatterDestroy(&matctx->scatter_id[i]));
1261: for (i=0;i<pep->nmat;i++) PetscCall(MatDestroy(&matctx->A[i]));
1262: if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
1263: for (i=0;i<matctx->subc->n;i++) PetscCall(VecScatterDestroy(&matctx->scatterp_id[i]));
1264: PetscCall(VecDestroy(&matctx->tp));
1265: PetscCall(VecDestroy(&matctx->tpg));
1266: PetscCall(BVDestroy(&matctx->W));
1267: }
1268: PetscCall(PetscFree3(matctx->A,matctx->scatter_id,matctx->scatterp_id));
1269: PetscCall(BVDestroy(&matctx->V));
1270: PetscCall(VecDestroy(&matctx->t));
1271: PetscCall(VecDestroy(&matctx->tg));
1272: PetscCall(VecDestroy(&matctx->Rv));
1273: PetscCall(VecDestroy(&matctx->Vi));
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1277: PetscErrorCode PEPNewtonRefinement_TOAR(PEP pep,PetscScalar sigma,PetscInt *maxits,PetscReal *tol,PetscInt k,PetscScalar *S,PetscInt lds)
1278: {
1279: PetscScalar *H,*work,*dH,*fH,*dVS;
1280: PetscInt ldh,i,j,its=1,nmat=pep->nmat,nsubc=pep->npart,rds;
1281: PetscBLASInt k_,ld_,*p;
1282: BV dV;
1283: PetscBool sinvert,flg;
1284: PEP_REFINE_EXPLICIT *matctx=NULL;
1285: Vec v;
1286: Mat M,P;
1287: PEP_REFINE_MATSHELL *ctx;
1289: PetscFunctionBegin;
1290: PetscCall(PetscLogEventBegin(PEP_Refine,pep,0,0,0));
1291: PetscCheck(k<=pep->n,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Multiple Refinement available only for invariant pairs of dimension smaller than n=%" PetscInt_FMT,pep->n);
1292: /* the input tolerance is not being taken into account (by the moment) */
1293: its = *maxits;
1294: PetscCall(PetscMalloc3(k*k,&dH,nmat*k*k,&fH,k,&work));
1295: PetscCall(DSGetLeadingDimension(pep->ds,&ldh));
1296: PetscCall(PetscMalloc1(2*k*k,&dVS));
1297: PetscCall(STGetTransform(pep->st,&flg));
1298: if (!flg && pep->st && pep->ops->backtransform) { /* BackTransform */
1299: PetscCall(PetscBLASIntCast(k,&k_));
1300: PetscCall(PetscBLASIntCast(ldh,&ld_));
1301: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert));
1302: if (sinvert) {
1303: PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
1304: PetscCall(PetscMalloc1(k,&p));
1305: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1306: PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&k_,&k_,H,&ld_,p,&info));
1307: PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&k_,H,&ld_,p,work,&k_,&info));
1308: PetscCall(PetscFPTrapPop());
1309: PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
1310: pep->ops->backtransform = NULL;
1311: }
1312: if (sigma!=0.0) {
1313: PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
1314: for (i=0;i<k;i++) H[i+ldh*i] += sigma;
1315: PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
1316: pep->ops->backtransform = NULL;
1317: }
1318: }
1319: if ((pep->scale==PEP_SCALE_BOTH || pep->scale==PEP_SCALE_SCALAR) && pep->sfactor!=1.0) {
1320: PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
1321: for (j=0;j<k;j++) {
1322: for (i=0;i<k;i++) H[i+j*ldh] *= pep->sfactor;
1323: }
1324: PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
1325: if (!flg) {
1326: /* Restore original values */
1327: for (i=0;i<pep->nmat;i++) {
1328: pep->pbc[pep->nmat+i] *= pep->sfactor;
1329: pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
1330: }
1331: }
1332: }
1333: if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr) {
1334: for (i=0;i<k;i++) {
1335: PetscCall(BVGetColumn(pep->V,i,&v));
1336: PetscCall(VecPointwiseMult(v,v,pep->Dr));
1337: PetscCall(BVRestoreColumn(pep->V,i,&v));
1338: }
1339: }
1340: PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
1342: PetscCall(NRefOrthogStep(pep,k,H,ldh,fH,S,lds));
1343: /* check if H is in Schur form */
1344: for (i=0;i<k-1;i++) {
1345: #if !PetscDefined(USE_COMPLEX)
1346: PetscCheck(H[i+1+i*ldh]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Iterative Refinement requires the complex Schur form of the projected matrix");
1347: #else
1348: PetscCheck(H[i+1+i*ldh]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Iterative Refinement requires an upper triangular projected matrix");
1349: #endif
1350: }
1351: PetscCheck(nsubc<=k,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Number of subcommunicators should not be larger than the invariant pair dimension");
1352: PetscCall(BVSetActiveColumns(pep->V,0,k));
1353: PetscCall(BVDuplicateResize(pep->V,k,&dV));
1354: if (pep->scheme!=PEP_REFINE_SCHEME_SCHUR) {
1355: PetscCall(PetscMalloc1(1,&matctx));
1356: if (nsubc>1) { /* splitting in subcommunicators */
1357: matctx->subc = pep->refinesubc;
1358: PetscCall(NRefSubcommSetup(pep,k,matctx,nsubc));
1359: } else matctx->subc=NULL;
1360: }
1362: /* Loop performing iterative refinements */
1363: for (i=0;i<its;i++) {
1364: /* Pre-compute the polynomial basis evaluated in H */
1365: PetscCall(PEPEvaluateBasisforMatrix(pep,nmat,k,H,ldh,fH));
1366: PetscCall(PEPNRefSetUp(pep,k,H,ldh,matctx,PetscNot(i)));
1367: /* Solve the linear system */
1368: PetscCall(PEPNRefForwardSubstitution(pep,k,S,lds,H,ldh,fH,dV,dVS,&rds,dH,k,pep->refineksp,matctx));
1369: /* Update X (=V*S) and H, and orthogonalize [X;X*fH1;...;XfH(deg-1)] */
1370: PetscCall(PEPNRefUpdateInvPair(pep,k,H,ldh,fH,dH,S,lds,dV,dVS,rds));
1371: }
1372: PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
1373: if (!flg && sinvert) PetscCall(PetscFree(p));
1374: PetscCall(PetscFree3(dH,fH,work));
1375: PetscCall(PetscFree(dVS));
1376: PetscCall(BVDestroy(&dV));
1377: switch (pep->scheme) {
1378: case PEP_REFINE_SCHEME_EXPLICIT:
1379: for (i=0;i<2;i++) PetscCall(MatDestroy(&matctx->E[i]));
1380: PetscCall(PetscFree4(matctx->idxp,matctx->idxg,matctx->map0,matctx->map1));
1381: PetscCall(VecDestroy(&matctx->tN));
1382: PetscCall(VecDestroy(&matctx->ttN));
1383: PetscCall(VecDestroy(&matctx->t1));
1384: if (nsubc>1) PetscCall(NRefSubcommDestroy(pep,matctx));
1385: else {
1386: PetscCall(VecDestroy(&matctx->vseq));
1387: PetscCall(VecScatterDestroy(&matctx->scatterctx));
1388: }
1389: PetscCall(PetscFree(matctx));
1390: PetscCall(KSPGetOperators(pep->refineksp,&M,NULL));
1391: PetscCall(MatDestroy(&M));
1392: break;
1393: case PEP_REFINE_SCHEME_MBE:
1394: PetscCall(BVDestroy(&matctx->W));
1395: PetscCall(BVDestroy(&matctx->Wt));
1396: PetscCall(BVDestroy(&matctx->M2));
1397: PetscCall(BVDestroy(&matctx->M3));
1398: PetscCall(MatDestroy(&matctx->M1));
1399: PetscCall(VecDestroy(&matctx->t));
1400: PetscCall(PetscFree5(matctx->M4,matctx->w,matctx->wt,matctx->d,matctx->dt));
1401: if (nsubc>1) PetscCall(NRefSubcommDestroy(pep,matctx));
1402: PetscCall(PetscFree(matctx));
1403: break;
1404: case PEP_REFINE_SCHEME_SCHUR:
1405: PetscCall(KSPGetOperators(pep->refineksp,&M,&P));
1406: PetscCall(MatShellGetContext(M,&ctx));
1407: PetscCall(PetscFree5(ctx->A,ctx->M4,ctx->pM4,ctx->work,ctx->fih));
1408: PetscCall(MatDestroy(&ctx->M1));
1409: PetscCall(BVDestroy(&ctx->M2));
1410: PetscCall(BVDestroy(&ctx->M3));
1411: PetscCall(BVDestroy(&ctx->W));
1412: PetscCall(VecDestroy(&ctx->t));
1413: PetscCall(PetscFree(ctx));
1414: PetscCall(MatDestroy(&M));
1415: PetscCall(MatDestroy(&P));
1416: break;
1417: }
1418: PetscCall(PetscLogEventEnd(PEP_Refine,pep,0,0,0));
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }