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