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 : SLEPc eigensolver: "davidson"
12 :
13 : Step: improve the eigenvectors X
14 : */
15 :
16 : #include "davidson.h"
17 : #include <slepcblaslapack.h>
18 :
19 : /**** JD update step (I - Kfg'/(g'Kf)) K(A - sB) (I - Kfg'/(g'Kf)) t = (I - Kfg'/(g'Kf))r ****/
20 :
21 : typedef struct {
22 : PetscInt size_X;
23 : KSP ksp; /* correction equation solver */
24 : Vec friends; /* reference vector for composite vectors */
25 : PetscScalar theta[4],thetai[2]; /* the shifts used in the correction eq. */
26 : PetscInt maxits; /* maximum number of iterations */
27 : PetscInt r_s,r_e; /* the selected eigenpairs to improve */
28 : PetscInt ksp_max_size; /* the ksp maximum subvectors size */
29 : PetscReal tol; /* the maximum solution tolerance */
30 : PetscReal lastTol; /* last tol for dynamic stopping criterion */
31 : PetscReal fix; /* tolerance for using the approx. eigenvalue */
32 : PetscBool dynamic; /* if the dynamic stopping criterion is applied */
33 : dvdDashboard *d; /* the current dvdDashboard reference */
34 : PC old_pc; /* old pc in ksp */
35 : BV KZ; /* KZ vecs for the projector KZ*inv(X'*KZ)*X' */
36 : BV U; /* new X vectors */
37 : PetscScalar *XKZ; /* X'*KZ */
38 : PetscScalar *iXKZ; /* inverse of XKZ */
39 : PetscInt ldXKZ; /* leading dimension of XKZ */
40 : PetscInt size_iXKZ; /* size of iXKZ */
41 : PetscInt ldiXKZ; /* leading dimension of iXKZ */
42 : PetscInt size_cX; /* last value of d->size_cX */
43 : PetscInt old_size_X; /* last number of improved vectors */
44 : PetscBLASInt *iXKZPivots; /* array of pivots */
45 : } dvdImprovex_jd;
46 :
47 : /*
48 : Compute (I - KZ*iXKZ*X')*V where,
49 : V, the vectors to apply the projector,
50 : cV, the number of vectors in V,
51 : */
52 13741 : static PetscErrorCode dvd_improvex_apply_proj(dvdDashboard *d,Vec *V,PetscInt cV)
53 : {
54 13741 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
55 13741 : PetscInt i,ldh,k,l;
56 13741 : PetscScalar *h;
57 13741 : PetscBLASInt cV_,n,info,ld;
58 : #if defined(PETSC_USE_COMPLEX)
59 : PetscInt j;
60 : #endif
61 :
62 13741 : PetscFunctionBegin;
63 13741 : PetscAssert(cV<=2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
64 :
65 : /* h <- X'*V */
66 13741 : PetscCall(PetscMalloc1(data->size_iXKZ*cV,&h));
67 13741 : ldh = data->size_iXKZ;
68 13741 : PetscCall(BVGetActiveColumns(data->U,&l,&k));
69 13741 : PetscAssert(ldh==k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
70 13741 : PetscCall(BVSetActiveColumns(data->U,0,k));
71 27703 : for (i=0;i<cV;i++) {
72 13962 : PetscCall(BVDotVec(data->U,V[i],&h[ldh*i]));
73 : #if defined(PETSC_USE_COMPLEX)
74 : for (j=0; j<k; j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
75 : #endif
76 : }
77 13741 : PetscCall(BVSetActiveColumns(data->U,l,k));
78 :
79 : /* h <- iXKZ\h */
80 13741 : PetscCall(PetscBLASIntCast(cV,&cV_));
81 13741 : PetscCall(PetscBLASIntCast(data->size_iXKZ,&n));
82 13741 : PetscCall(PetscBLASIntCast(data->ldiXKZ,&ld));
83 13741 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
84 13741 : PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
85 13741 : PetscCall(PetscFPTrapPop());
86 13741 : SlepcCheckLapackInfo("getrs",info);
87 :
88 : /* V <- V - KZ*h */
89 13741 : PetscCall(BVSetActiveColumns(data->KZ,0,k));
90 27703 : for (i=0;i<cV;i++) PetscCall(BVMultVec(data->KZ,-1.0,1.0,V[i],&h[ldh*i]));
91 13741 : PetscCall(BVSetActiveColumns(data->KZ,l,k));
92 13741 : PetscCall(PetscFree(h));
93 13741 : PetscFunctionReturn(PETSC_SUCCESS);
94 : }
95 :
96 : /*
97 : Compute (I - X*iXKZ*KZ')*V where,
98 : V, the vectors to apply the projector,
99 : cV, the number of vectors in V,
100 : */
101 91 : static PetscErrorCode dvd_improvex_applytrans_proj(dvdDashboard *d,Vec *V,PetscInt cV)
102 : {
103 91 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
104 91 : PetscInt i,ldh,k,l;
105 91 : PetscScalar *h;
106 91 : PetscBLASInt cV_, n, info, ld;
107 : #if defined(PETSC_USE_COMPLEX)
108 : PetscInt j;
109 : #endif
110 :
111 91 : PetscFunctionBegin;
112 91 : PetscAssert(cV<=2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
113 :
114 : /* h <- KZ'*V */
115 91 : PetscCall(PetscMalloc1(data->size_iXKZ*cV,&h));
116 91 : ldh = data->size_iXKZ;
117 91 : PetscCall(BVGetActiveColumns(data->U,&l,&k));
118 91 : PetscAssert(ldh==k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
119 91 : PetscCall(BVSetActiveColumns(data->KZ,0,k));
120 182 : for (i=0;i<cV;i++) {
121 91 : PetscCall(BVDotVec(data->KZ,V[i],&h[ldh*i]));
122 : #if defined(PETSC_USE_COMPLEX)
123 : for (j=0;j<k;j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
124 : #endif
125 : }
126 91 : PetscCall(BVSetActiveColumns(data->KZ,l,k));
127 :
128 : /* h <- iXKZ\h */
129 91 : PetscCall(PetscBLASIntCast(cV,&cV_));
130 91 : PetscCall(PetscBLASIntCast(data->size_iXKZ,&n));
131 91 : PetscCall(PetscBLASIntCast(data->ldiXKZ,&ld));
132 91 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
133 91 : PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
134 91 : PetscCall(PetscFPTrapPop());
135 91 : SlepcCheckLapackInfo("getrs",info);
136 :
137 : /* V <- V - U*h */
138 91 : PetscCall(BVSetActiveColumns(data->U,0,k));
139 182 : for (i=0;i<cV;i++) PetscCall(BVMultVec(data->U,-1.0,1.0,V[i],&h[ldh*i]));
140 91 : PetscCall(BVSetActiveColumns(data->U,l,k));
141 91 : PetscCall(PetscFree(h));
142 91 : PetscFunctionReturn(PETSC_SUCCESS);
143 : }
144 :
145 78 : static PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)
146 : {
147 78 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
148 :
149 78 : PetscFunctionBegin;
150 78 : PetscCall(VecDestroy(&data->friends));
151 :
152 : /* Restore the pc of ksp */
153 78 : if (data->old_pc) {
154 18 : PetscCall(KSPSetPC(data->ksp, data->old_pc));
155 18 : PetscCall(PCDestroy(&data->old_pc));
156 : }
157 78 : PetscFunctionReturn(PETSC_SUCCESS);
158 : }
159 :
160 78 : static PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)
161 : {
162 78 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
163 :
164 78 : PetscFunctionBegin;
165 : /* Free local data and objects */
166 78 : PetscCall(PetscFree(data->XKZ));
167 78 : PetscCall(PetscFree(data->iXKZ));
168 78 : PetscCall(PetscFree(data->iXKZPivots));
169 78 : PetscCall(BVDestroy(&data->KZ));
170 78 : PetscCall(BVDestroy(&data->U));
171 78 : PetscCall(PetscFree(data));
172 78 : PetscFunctionReturn(PETSC_SUCCESS);
173 : }
174 :
175 : /*
176 : y <- theta[1]A*x - theta[0]*B*x
177 : auxV, two auxiliary vectors
178 : */
179 18526 : static inline PetscErrorCode dvd_aux_matmult(dvdImprovex_jd *data,const Vec *x,const Vec *y)
180 : {
181 18526 : PetscInt n,i;
182 18526 : const Vec *Bx;
183 18526 : Vec *auxV;
184 :
185 18526 : PetscFunctionBegin;
186 18526 : n = data->r_e - data->r_s;
187 37224 : for (i=0;i<n;i++) PetscCall(MatMult(data->d->A,x[i],y[i]));
188 :
189 18526 : PetscCall(SlepcVecPoolGetVecs(data->d->auxV,2,&auxV));
190 37052 : for (i=0;i<n;i++) {
191 : #if !defined(PETSC_USE_COMPLEX)
192 18526 : if (PetscUnlikely(data->d->eigi[data->r_s+i] != 0.0)) {
193 172 : if (data->d->B) {
194 0 : PetscCall(MatMult(data->d->B,x[i],auxV[0]));
195 0 : PetscCall(MatMult(data->d->B,x[i+1],auxV[1]));
196 0 : Bx = auxV;
197 172 : } else Bx = &x[i];
198 :
199 : /* y_i <- [ t_2i+1*A*x_i - t_2i*Bx_i + ti_i*Bx_i+1;
200 : y_i+1 t_2i+1*A*x_i+1 - ti_i*Bx_i - t_2i*Bx_i+1 ] */
201 172 : PetscCall(VecAXPBYPCZ(y[i],-data->theta[2*i],data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]));
202 172 : PetscCall(VecAXPBYPCZ(y[i+1],-data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]));
203 172 : i++;
204 : } else
205 : #endif
206 : {
207 18354 : if (data->d->B) {
208 3400 : PetscCall(MatMult(data->d->B,x[i],auxV[0]));
209 3400 : Bx = auxV;
210 14954 : } else Bx = &x[i];
211 18526 : PetscCall(VecAXPBY(y[i],-data->theta[i*2],data->theta[i*2+1],Bx[0]));
212 : }
213 : }
214 18526 : PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV));
215 18526 : PetscFunctionReturn(PETSC_SUCCESS);
216 : }
217 :
218 : /*
219 : y <- theta[1]'*A'*x - theta[0]'*B'*x
220 : */
221 75 : static inline PetscErrorCode dvd_aux_matmulttrans(dvdImprovex_jd *data,const Vec *x,const Vec *y)
222 : {
223 75 : PetscInt n,i;
224 75 : const Vec *Bx;
225 75 : Vec *auxV;
226 :
227 75 : PetscFunctionBegin;
228 75 : n = data->r_e - data->r_s;
229 150 : for (i=0;i<n;i++) PetscCall(MatMultTranspose(data->d->A,x[i],y[i]));
230 :
231 75 : PetscCall(SlepcVecPoolGetVecs(data->d->auxV,2,&auxV));
232 150 : for (i=0;i<n;i++) {
233 : #if !defined(PETSC_USE_COMPLEX)
234 75 : if (data->d->eigi[data->r_s+i] != 0.0) {
235 0 : if (data->d->B) {
236 0 : PetscCall(MatMultTranspose(data->d->B,x[i],auxV[0]));
237 0 : PetscCall(MatMultTranspose(data->d->B,x[i+1],auxV[1]));
238 0 : Bx = auxV;
239 0 : } else Bx = &x[i];
240 :
241 : /* y_i <- [ t_2i+1*A*x_i - t_2i*Bx_i - ti_i*Bx_i+1;
242 : y_i+1 t_2i+1*A*x_i+1 + ti_i*Bx_i - t_2i*Bx_i+1 ] */
243 0 : PetscCall(VecAXPBYPCZ(y[i],-data->theta[2*i],-data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]));
244 0 : PetscCall(VecAXPBYPCZ(y[i+1],data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]));
245 0 : i++;
246 : } else
247 : #endif
248 : {
249 75 : if (data->d->B) {
250 0 : PetscCall(MatMultTranspose(data->d->B,x[i],auxV[0]));
251 0 : Bx = auxV;
252 75 : } else Bx = &x[i];
253 75 : PetscCall(VecAXPBY(y[i],PetscConj(-data->theta[i*2]),PetscConj(data->theta[i*2+1]),Bx[0]));
254 : }
255 : }
256 75 : PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV));
257 75 : PetscFunctionReturn(PETSC_SUCCESS);
258 : }
259 :
260 9266 : static PetscErrorCode PCApplyBA_dvd(PC pc,PCSide side,Vec in,Vec out,Vec w)
261 : {
262 9266 : dvdImprovex_jd *data;
263 9266 : PetscInt n,i;
264 9266 : const Vec *inx,*outx,*wx;
265 9266 : Vec *auxV;
266 9266 : Mat A;
267 :
268 9266 : PetscFunctionBegin;
269 9266 : PetscCall(PCGetOperators(pc,&A,NULL));
270 9266 : PetscCall(MatShellGetContext(A,&data));
271 9266 : PetscCall(VecCompGetSubVecs(in,NULL,&inx));
272 9266 : PetscCall(VecCompGetSubVecs(out,NULL,&outx));
273 9266 : PetscCall(VecCompGetSubVecs(w,NULL,&wx));
274 9266 : n = data->r_e - data->r_s;
275 9266 : PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
276 9266 : switch (side) {
277 9266 : case PC_LEFT:
278 : /* aux <- theta[1]A*in - theta[0]*B*in */
279 9266 : PetscCall(dvd_aux_matmult(data,inx,auxV));
280 :
281 : /* out <- K * aux */
282 18704 : for (i=0;i<n;i++) PetscCall(data->d->improvex_precond(data->d,data->r_s+i,auxV[i],outx[i]));
283 : break;
284 : case PC_RIGHT:
285 : /* aux <- K * in */
286 0 : for (i=0;i<n;i++) PetscCall(data->d->improvex_precond(data->d,data->r_s+i,inx[i],auxV[i]));
287 :
288 : /* out <- theta[1]A*auxV - theta[0]*B*auxV */
289 0 : PetscCall(dvd_aux_matmult(data,auxV,outx));
290 : break;
291 : case PC_SYMMETRIC:
292 : /* aux <- K^{1/2} * in */
293 0 : for (i=0;i<n;i++) PetscCall(PCApplySymmetricRight(data->old_pc,inx[i],auxV[i]));
294 :
295 : /* wx <- theta[1]A*auxV - theta[0]*B*auxV */
296 0 : PetscCall(dvd_aux_matmult(data,auxV,wx));
297 :
298 : /* aux <- K^{1/2} * in */
299 0 : for (i=0;i<n;i++) PetscCall(PCApplySymmetricLeft(data->old_pc,wx[i],outx[i]));
300 : break;
301 0 : default:
302 0 : SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported KSP side");
303 : }
304 : /* out <- out - v*(u'*out) */
305 9266 : PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
306 9266 : PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
307 9266 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 :
310 381 : static PetscErrorCode PCApply_dvd(PC pc,Vec in,Vec out)
311 : {
312 381 : dvdImprovex_jd *data;
313 381 : PetscInt n,i;
314 381 : const Vec *inx, *outx;
315 381 : Mat A;
316 :
317 381 : PetscFunctionBegin;
318 381 : PetscCall(PCGetOperators(pc,&A,NULL));
319 381 : PetscCall(MatShellGetContext(A,&data));
320 381 : PetscCall(VecCompGetSubVecs(in,NULL,&inx));
321 381 : PetscCall(VecCompGetSubVecs(out,NULL,&outx));
322 381 : n = data->r_e - data->r_s;
323 : /* out <- K * in */
324 765 : for (i=0;i<n;i++) PetscCall(data->d->improvex_precond(data->d,data->r_s+i,inx[i],outx[i]));
325 : /* out <- out - v*(u'*out) */
326 381 : PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
327 381 : PetscFunctionReturn(PETSC_SUCCESS);
328 : }
329 :
330 91 : static PetscErrorCode PCApplyTranspose_dvd(PC pc,Vec in,Vec out)
331 : {
332 91 : dvdImprovex_jd *data;
333 91 : PetscInt n,i;
334 91 : const Vec *inx, *outx;
335 91 : Vec *auxV;
336 91 : Mat A;
337 :
338 91 : PetscFunctionBegin;
339 91 : PetscCall(PCGetOperators(pc,&A,NULL));
340 91 : PetscCall(MatShellGetContext(A,&data));
341 91 : PetscCall(VecCompGetSubVecs(in,NULL,&inx));
342 91 : PetscCall(VecCompGetSubVecs(out,NULL,&outx));
343 91 : n = data->r_e - data->r_s;
344 91 : PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
345 : /* auxV <- in */
346 182 : for (i=0;i<n;i++) PetscCall(VecCopy(inx[i],auxV[i]));
347 : /* auxV <- auxV - u*(v'*auxV) */
348 91 : PetscCall(dvd_improvex_applytrans_proj(data->d,auxV,n));
349 : /* out <- K' * aux */
350 182 : for (i=0;i<n;i++) PetscCall(PCApplyTranspose(data->old_pc,auxV[i],outx[i]));
351 91 : PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
352 91 : PetscFunctionReturn(PETSC_SUCCESS);
353 : }
354 :
355 9260 : static PetscErrorCode MatMult_dvd_jd(Mat A,Vec in,Vec out)
356 : {
357 9260 : dvdImprovex_jd *data;
358 9260 : PetscInt n;
359 9260 : const Vec *inx, *outx;
360 9260 : PCSide side;
361 :
362 9260 : PetscFunctionBegin;
363 9260 : PetscCall(MatShellGetContext(A,&data));
364 9260 : PetscCall(VecCompGetSubVecs(in,NULL,&inx));
365 9260 : PetscCall(VecCompGetSubVecs(out,NULL,&outx));
366 9260 : n = data->r_e - data->r_s;
367 : /* out <- theta[1]A*in - theta[0]*B*in */
368 9260 : PetscCall(dvd_aux_matmult(data,inx,outx));
369 9260 : PetscCall(KSPGetPCSide(data->ksp,&side));
370 9260 : if (side == PC_RIGHT) {
371 : /* out <- out - v*(u'*out) */
372 0 : PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
373 : }
374 9260 : PetscFunctionReturn(PETSC_SUCCESS);
375 : }
376 :
377 75 : static PetscErrorCode MatMultTranspose_dvd_jd(Mat A,Vec in,Vec out)
378 : {
379 75 : dvdImprovex_jd *data;
380 75 : PetscInt n,i;
381 75 : const Vec *inx,*outx,*r;
382 75 : Vec *auxV;
383 75 : PCSide side;
384 :
385 75 : PetscFunctionBegin;
386 75 : PetscCall(MatShellGetContext(A,&data));
387 75 : PetscCall(VecCompGetSubVecs(in,NULL,&inx));
388 75 : PetscCall(VecCompGetSubVecs(out,NULL,&outx));
389 75 : n = data->r_e - data->r_s;
390 75 : PetscCall(KSPGetPCSide(data->ksp,&side));
391 75 : if (side == PC_RIGHT) {
392 : /* auxV <- in */
393 0 : PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
394 0 : for (i=0;i<n;i++) PetscCall(VecCopy(inx[i],auxV[i]));
395 : /* auxV <- auxV - v*(u'*auxV) */
396 0 : PetscCall(dvd_improvex_applytrans_proj(data->d,auxV,n));
397 0 : r = auxV;
398 75 : } else r = inx;
399 : /* out <- theta[1]A*r - theta[0]*B*r */
400 75 : PetscCall(dvd_aux_matmulttrans(data,r,outx));
401 75 : if (side == PC_RIGHT) PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
402 75 : PetscFunctionReturn(PETSC_SUCCESS);
403 : }
404 :
405 24 : static PetscErrorCode MatCreateVecs_dvd_jd(Mat A,Vec *right,Vec *left)
406 : {
407 24 : Vec *r,*l;
408 24 : dvdImprovex_jd *data;
409 24 : PetscInt n,i;
410 :
411 24 : PetscFunctionBegin;
412 24 : PetscCall(MatShellGetContext(A,&data));
413 24 : n = data->ksp_max_size;
414 24 : if (right) PetscCall(PetscMalloc1(n,&r));
415 24 : if (left) PetscCall(PetscMalloc1(n,&l));
416 56 : for (i=0;i<n;i++) PetscCall(MatCreateVecs(data->d->A,right?&r[i]:NULL,left?&l[i]:NULL));
417 24 : if (right) {
418 24 : PetscCall(VecCreateCompWithVecs(r,n,data->friends,right));
419 56 : for (i=0;i<n;i++) PetscCall(VecDestroy(&r[i]));
420 : }
421 24 : if (left) {
422 0 : PetscCall(VecCreateCompWithVecs(l,n,data->friends,left));
423 0 : for (i=0;i<n;i++) PetscCall(VecDestroy(&l[i]));
424 : }
425 :
426 24 : if (right) PetscCall(PetscFree(r));
427 24 : if (left) PetscCall(PetscFree(l));
428 24 : PetscFunctionReturn(PETSC_SUCCESS);
429 : }
430 :
431 78 : static PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)
432 : {
433 78 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
434 78 : PetscInt rA, cA, rlA, clA;
435 78 : Mat A;
436 78 : PetscBool t;
437 78 : PC pc;
438 78 : Vec v0[2];
439 :
440 78 : PetscFunctionBegin;
441 78 : data->size_cX = data->old_size_X = 0;
442 78 : data->lastTol = data->dynamic?0.5:0.0;
443 :
444 : /* Setup the ksp */
445 78 : if (data->ksp) {
446 : /* Create the reference vector */
447 31 : PetscCall(BVGetColumn(d->eps->V,0,&v0[0]));
448 31 : v0[1] = v0[0];
449 31 : PetscCall(VecCreateCompWithVecs(v0,data->ksp_max_size,NULL,&data->friends));
450 31 : PetscCall(BVRestoreColumn(d->eps->V,0,&v0[0]));
451 :
452 : /* Save the current pc and set a PCNONE */
453 31 : PetscCall(KSPGetPC(data->ksp, &data->old_pc));
454 31 : PetscCall(PetscObjectTypeCompare((PetscObject)data->old_pc,PCNONE,&t));
455 31 : data->lastTol = 0.5;
456 31 : if (t) data->old_pc = NULL;
457 : else {
458 18 : PetscCall(PetscObjectReference((PetscObject)data->old_pc));
459 18 : PetscCall(PCCreate(PetscObjectComm((PetscObject)d->eps),&pc));
460 18 : PetscCall(PCSetType(pc,PCSHELL));
461 18 : PetscCall(PCSetOperators(pc,d->A,d->A));
462 18 : PetscCall(PCSetReusePreconditioner(pc,PETSC_TRUE));
463 18 : PetscCall(PCShellSetApply(pc,PCApply_dvd));
464 18 : PetscCall(PCShellSetApplyBA(pc,PCApplyBA_dvd));
465 18 : PetscCall(PCShellSetApplyTranspose(pc,PCApplyTranspose_dvd));
466 18 : PetscCall(KSPSetPC(data->ksp,pc));
467 18 : PetscCall(PCDestroy(&pc));
468 : }
469 :
470 : /* Create the (I-v*u')*K*(A-s*B) matrix */
471 31 : PetscCall(MatGetSize(d->A,&rA,&cA));
472 31 : PetscCall(MatGetLocalSize(d->A,&rlA,&clA));
473 31 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)d->A),rlA*data->ksp_max_size,clA*data->ksp_max_size,rA*data->ksp_max_size,cA*data->ksp_max_size,data,&A));
474 31 : PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_dvd_jd));
475 31 : PetscCall(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_dvd_jd));
476 31 : PetscCall(MatShellSetOperation(A,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_dvd_jd));
477 :
478 : /* Try to avoid KSPReset */
479 31 : PetscCall(KSPGetOperatorsSet(data->ksp,&t,NULL));
480 31 : if (t) {
481 24 : Mat M;
482 24 : PetscInt rM;
483 24 : PetscCall(KSPGetOperators(data->ksp,&M,NULL));
484 24 : PetscCall(MatGetSize(M,&rM,NULL));
485 24 : if (rM != rA*data->ksp_max_size) PetscCall(KSPReset(data->ksp));
486 : }
487 31 : PetscCall(EPS_KSPSetOperators(data->ksp,A,A));
488 31 : PetscCall(KSPSetReusePreconditioner(data->ksp,PETSC_TRUE));
489 31 : PetscCall(KSPSetUp(data->ksp));
490 31 : PetscCall(MatDestroy(&A));
491 : } else {
492 47 : data->old_pc = NULL;
493 47 : data->friends = NULL;
494 : }
495 78 : PetscCall(BVSetActiveColumns(data->KZ,0,0));
496 78 : PetscCall(BVSetActiveColumns(data->U,0,0));
497 78 : PetscFunctionReturn(PETSC_SUCCESS);
498 : }
499 :
500 : /*
501 : Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
502 : kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
503 : where
504 : pX,pY, the right and left eigenvectors of the projected system
505 : ld, the leading dimension of pX and pY
506 : */
507 4826 : static PetscErrorCode dvd_improvex_jd_proj_cuv(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *kr,PetscScalar *theta,PetscScalar *thetai,PetscScalar *pX,PetscScalar *pY,PetscInt ld)
508 : {
509 4826 : PetscInt n=i_e-i_s,size_KZ,V_new,rm,i,lv,kv,lKZ,kKZ;
510 4826 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
511 4826 : const PetscScalar *array;
512 4826 : Mat M;
513 4826 : Vec u[2],v[2];
514 4826 : PetscBLASInt s,ldXKZ,info;
515 :
516 4826 : PetscFunctionBegin;
517 : /* Check consistency */
518 4826 : PetscCall(BVGetActiveColumns(d->eps->V,&lv,&kv));
519 4826 : V_new = lv - data->size_cX;
520 4826 : PetscAssert(V_new<=data->old_size_X,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
521 4826 : data->old_size_X = n;
522 4826 : data->size_cX = lv;
523 :
524 : /* KZ <- KZ(rm:rm+max_cX-1) */
525 4826 : PetscCall(BVGetActiveColumns(data->KZ,&lKZ,&kKZ));
526 4826 : rm = PetscMax(V_new+lKZ,0);
527 4826 : if (rm > 0) {
528 175 : for (i=0;i<lKZ;i++) {
529 0 : PetscCall(BVCopyColumn(data->KZ,i+rm,i));
530 0 : PetscCall(BVCopyColumn(data->U,i+rm,i));
531 : }
532 : }
533 :
534 : /* XKZ <- XKZ(rm:rm+max_cX-1,rm:rm+max_cX-1) */
535 4826 : if (rm > 0) {
536 175 : for (i=0;i<lKZ;i++) PetscCall(PetscArraycpy(&data->XKZ[i*data->ldXKZ+i],&data->XKZ[(i+rm)*data->ldXKZ+i+rm],lKZ));
537 : }
538 4826 : lKZ = PetscMin(0,lKZ+V_new);
539 4826 : PetscCall(BVSetActiveColumns(data->KZ,lKZ,lKZ+n));
540 4826 : PetscCall(BVSetActiveColumns(data->U,lKZ,lKZ+n));
541 :
542 : /* Compute X, KZ and KR */
543 4826 : PetscCall(BVGetColumn(data->U,lKZ,u));
544 4826 : if (n>1) PetscCall(BVGetColumn(data->U,lKZ+1,&u[1]));
545 4826 : PetscCall(BVGetColumn(data->KZ,lKZ,v));
546 4826 : if (n>1) PetscCall(BVGetColumn(data->KZ,lKZ+1,&v[1]));
547 4826 : PetscCall(d->improvex_jd_proj_uv(d,i_s,i_e,u,v,kr,theta,thetai,pX,pY,ld));
548 4826 : PetscCall(BVRestoreColumn(data->U,lKZ,u));
549 4826 : if (n>1) PetscCall(BVRestoreColumn(data->U,lKZ+1,&u[1]));
550 4826 : PetscCall(BVRestoreColumn(data->KZ,lKZ,v));
551 4826 : if (n>1) PetscCall(BVRestoreColumn(data->KZ,lKZ+1,&v[1]));
552 :
553 : /* XKZ <- U'*KZ */
554 4826 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,lKZ+n,lKZ+n,NULL,&M));
555 4826 : PetscCall(BVMatProject(data->KZ,NULL,data->U,M));
556 4826 : PetscCall(MatDenseGetArrayRead(M,&array));
557 9703 : for (i=lKZ;i<lKZ+n;i++) { /* upper part */
558 4877 : PetscCall(PetscArraycpy(&data->XKZ[data->ldXKZ*i],&array[i*(lKZ+n)],lKZ));
559 : }
560 9703 : for (i=0;i<lKZ+n;i++) { /* lower part */
561 4877 : PetscCall(PetscArraycpy(&data->XKZ[data->ldXKZ*i+lKZ],&array[i*(lKZ+n)+lKZ],n));
562 : }
563 4826 : PetscCall(MatDenseRestoreArrayRead(M,&array));
564 4826 : PetscCall(MatDestroy(&M));
565 :
566 : /* iXKZ <- inv(XKZ) */
567 4826 : size_KZ = lKZ+n;
568 4826 : PetscCall(PetscBLASIntCast(lKZ+n,&s));
569 4826 : data->ldiXKZ = data->size_iXKZ = size_KZ;
570 9703 : for (i=0;i<size_KZ;i++) PetscCall(PetscArraycpy(&data->iXKZ[data->ldiXKZ*i],&data->XKZ[data->ldXKZ*i],size_KZ));
571 4826 : PetscCall(PetscBLASIntCast(data->ldiXKZ,&ldXKZ));
572 4826 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
573 4826 : PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&s,&s,data->iXKZ,&ldXKZ,data->iXKZPivots,&info));
574 4826 : PetscCall(PetscFPTrapPop());
575 4826 : SlepcCheckLapackInfo("getrf",info);
576 4826 : PetscFunctionReturn(PETSC_SUCCESS);
577 : }
578 :
579 4446 : static PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
580 : {
581 4446 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
582 4446 : PetscInt i,j,n,maxits,maxits0,lits,s,ld,k,max_size_D,lV,kV;
583 4446 : PetscScalar *pX,*pY;
584 4446 : PetscReal tol,tol0;
585 4446 : Vec *kr,kr_comp,D_comp,D[2],kr0[2];
586 4446 : PetscBool odd_situation = PETSC_FALSE;
587 :
588 4446 : PetscFunctionBegin;
589 4446 : PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
590 4446 : max_size_D = d->eps->ncv-kV;
591 : /* Quick exit */
592 4446 : if ((max_size_D == 0) || r_e-r_s <= 0) {
593 0 : *size_D = 0;
594 0 : PetscFunctionReturn(PETSC_SUCCESS);
595 : }
596 :
597 4446 : n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);
598 4446 : PetscAssert(n>0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"n == 0");
599 4446 : PetscAssert(data->size_X>=r_e-r_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"size_X < r_e-r_s");
600 :
601 4446 : PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
602 :
603 : /* Restart lastTol if a new pair converged */
604 4446 : if (data->dynamic && data->size_cX < lV)
605 1 : data->lastTol = 0.5;
606 :
607 9019 : for (i=0,s=0;i<n;i+=s) {
608 : /* If the selected eigenvalue is complex, but the arithmetic is real... */
609 : #if !defined(PETSC_USE_COMPLEX)
610 4826 : if (d->eigi[r_s+i] != 0.0) {
611 51 : if (i+2 <= max_size_D) s=2;
612 : else break;
613 : } else
614 : #endif
615 : s=1;
616 :
617 4826 : data->r_s = r_s+i;
618 4826 : data->r_e = r_s+i+s;
619 4826 : PetscCall(SlepcVecPoolGetVecs(d->auxV,s,&kr));
620 :
621 : /* Compute theta, maximum iterations and tolerance */
622 : maxits = 0;
623 : tol = 1;
624 9703 : for (j=0;j<s;j++) {
625 4877 : PetscCall(d->improvex_jd_lit(d,r_s+i+j,&data->theta[2*j],&data->thetai[j],&maxits0,&tol0));
626 4877 : maxits += maxits0;
627 4877 : tol *= tol0;
628 : }
629 4826 : maxits/= s;
630 4826 : tol = data->dynamic?data->lastTol:PetscExpReal(PetscLogReal(tol)/s);
631 :
632 : /* Compute u, v and kr */
633 4826 : k = r_s+i;
634 4826 : PetscCall(DSVectors(d->eps->ds,DS_MAT_X,&k,NULL));
635 4826 : k = r_s+i;
636 4826 : PetscCall(DSVectors(d->eps->ds,DS_MAT_Y,&k,NULL));
637 4826 : PetscCall(DSGetArray(d->eps->ds,DS_MAT_X,&pX));
638 4826 : PetscCall(DSGetArray(d->eps->ds,DS_MAT_Y,&pY));
639 4826 : PetscCall(dvd_improvex_jd_proj_cuv(d,r_s+i,r_s+i+s,kr,data->theta,data->thetai,pX,pY,ld));
640 4826 : PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_X,&pX));
641 4826 : PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_Y,&pY));
642 :
643 : /* Check if the first eigenpairs are converged */
644 4826 : if (i == 0) {
645 4446 : PetscInt oldnpreconv = d->npreconv;
646 4446 : PetscCall(d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv));
647 4446 : if (d->npreconv > oldnpreconv) break;
648 : }
649 :
650 : /* Test the odd situation of solving Ax=b with A=I */
651 : #if !defined(PETSC_USE_COMPLEX)
652 4573 : odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && data->thetai[0] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
653 : #else
654 : odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
655 : #endif
656 : /* If JD */
657 4573 : if (data->ksp && !odd_situation) {
658 : /* kr <- -kr */
659 961 : for (j=0;j<s;j++) PetscCall(VecScale(kr[j],-1.0));
660 :
661 : /* Compose kr and D */
662 479 : kr0[0] = kr[0];
663 479 : kr0[1] = (s==2 ? kr[1] : NULL);
664 479 : PetscCall(VecCreateCompWithVecs(kr0,data->ksp_max_size,data->friends,&kr_comp));
665 479 : PetscCall(BVGetColumn(d->eps->V,kV+i,&D[0]));
666 479 : if (s==2) PetscCall(BVGetColumn(d->eps->V,kV+i+1,&D[1]));
667 476 : else D[1] = NULL;
668 479 : PetscCall(VecCreateCompWithVecs(D,data->ksp_max_size,data->friends,&D_comp));
669 479 : PetscCall(VecCompSetSubVecs(data->friends,s,NULL));
670 :
671 : /* Solve the correction equation */
672 479 : PetscCall(KSPSetTolerances(data->ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,maxits));
673 479 : PetscCall(KSPSolve(data->ksp,kr_comp,D_comp));
674 479 : PetscCall(KSPGetIterationNumber(data->ksp,&lits));
675 :
676 : /* Destroy the composed ks and D */
677 479 : PetscCall(VecDestroy(&kr_comp));
678 479 : PetscCall(VecDestroy(&D_comp));
679 479 : PetscCall(BVRestoreColumn(d->eps->V,kV+i,&D[0]));
680 479 : if (s==2) PetscCall(BVRestoreColumn(d->eps->V,kV+i+1,&D[1]));
681 :
682 : /* If GD */
683 : } else {
684 4094 : PetscCall(BVGetColumn(d->eps->V,kV+i,&D[0]));
685 4094 : if (s==2) PetscCall(BVGetColumn(d->eps->V,kV+i+1,&D[1]));
686 8234 : for (j=0;j<s;j++) PetscCall(d->improvex_precond(d,r_s+i+j,kr[j],D[j]));
687 4094 : PetscCall(dvd_improvex_apply_proj(d,D,s));
688 4094 : PetscCall(BVRestoreColumn(d->eps->V,kV+i,&D[0]));
689 4094 : if (s==2) PetscCall(BVRestoreColumn(d->eps->V,kV+i+1,&D[1]));
690 : }
691 : /* Prevent that short vectors are discarded in the orthogonalization */
692 4573 : if (i == 0 && d->eps->errest[d->nconv+r_s] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s] < PETSC_MAX_REAL) {
693 8434 : for (j=0;j<s;j++) PetscCall(BVScaleColumn(d->eps->V,kV+i+j,1.0/d->eps->errest[d->nconv+r_s]));
694 : }
695 4573 : PetscCall(SlepcVecPoolRestoreVecs(d->auxV,s,&kr));
696 : }
697 4446 : *size_D = i;
698 4467 : if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);
699 4446 : PetscFunctionReturn(PETSC_SUCCESS);
700 : }
701 :
702 234 : PetscErrorCode dvd_improvex_jd(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs,PetscBool dynamic)
703 : {
704 234 : dvdImprovex_jd *data;
705 234 : PetscBool useGD;
706 234 : PC pc;
707 234 : PetscInt size_P;
708 :
709 234 : PetscFunctionBegin;
710 : /* Setting configuration constrains */
711 234 : PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&useGD));
712 :
713 : /* If the arithmetic is real and the problem is not Hermitian, then
714 : the block size is incremented in one */
715 : #if !defined(PETSC_USE_COMPLEX)
716 234 : if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) {
717 66 : max_bs++;
718 66 : b->max_size_P = PetscMax(b->max_size_P,2);
719 : } else
720 : #endif
721 : {
722 168 : b->max_size_P = PetscMax(b->max_size_P,1);
723 : }
724 234 : b->max_size_X = PetscMax(b->max_size_X,max_bs);
725 234 : size_P = b->max_size_P;
726 :
727 : /* Setup the preconditioner */
728 234 : if (ksp) {
729 234 : PetscCall(KSPGetPC(ksp,&pc));
730 234 : PetscCall(dvd_static_precond_PC(d,b,pc));
731 0 : } else PetscCall(dvd_static_precond_PC(d,b,NULL));
732 :
733 : /* Setup the step */
734 234 : if (b->state >= DVD_STATE_CONF) {
735 78 : PetscCall(PetscNew(&data));
736 78 : data->dynamic = dynamic;
737 78 : PetscCall(PetscMalloc1(size_P*size_P,&data->XKZ));
738 78 : PetscCall(PetscMalloc1(size_P*size_P,&data->iXKZ));
739 78 : PetscCall(PetscMalloc1(size_P,&data->iXKZPivots));
740 78 : data->ldXKZ = size_P;
741 78 : data->size_X = b->max_size_X;
742 78 : d->improveX_data = data;
743 78 : data->ksp = useGD? NULL: ksp;
744 78 : data->d = d;
745 78 : d->improveX = dvd_improvex_jd_gen;
746 : #if !defined(PETSC_USE_COMPLEX)
747 78 : if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) data->ksp_max_size = 2;
748 : else
749 : #endif
750 56 : data->ksp_max_size = 1;
751 : /* Create various vector basis */
752 78 : PetscCall(BVDuplicateResize(d->eps->V,size_P,&data->KZ));
753 78 : PetscCall(BVSetMatrix(data->KZ,NULL,PETSC_FALSE));
754 78 : PetscCall(BVDuplicate(data->KZ,&data->U));
755 :
756 78 : PetscCall(EPSDavidsonFLAdd(&d->startList,dvd_improvex_jd_start));
757 78 : PetscCall(EPSDavidsonFLAdd(&d->endList,dvd_improvex_jd_end));
758 78 : PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_jd_d));
759 : }
760 234 : PetscFunctionReturn(PETSC_SUCCESS);
761 : }
762 :
763 : #if !defined(PETSC_USE_COMPLEX)
764 51 : static inline PetscErrorCode dvd_complex_rayleigh_quotient(Vec ur,Vec ui,Vec Axr,Vec Axi,Vec Bxr,Vec Bxi,PetscScalar *eigr,PetscScalar *eigi)
765 : {
766 51 : PetscScalar rAr,iAr,rAi,iAi,rBr,iBr,rBi,iBi,b0,b2,b4,b6,b7;
767 :
768 51 : PetscFunctionBegin;
769 : /* eigr = [(rAr+iAi)*(rBr+iBi) + (rAi-iAr)*(rBi-iBr)]/k
770 : eigi = [(rAi-iAr)*(rBr+iBi) - (rAr+iAi)*(rBi-iBr)]/k
771 : k = (rBr+iBi)*(rBr+iBi) + (rBi-iBr)*(rBi-iBr) */
772 51 : PetscCall(VecDotBegin(Axr,ur,&rAr)); /* r*A*r */
773 51 : PetscCall(VecDotBegin(Axr,ui,&iAr)); /* i*A*r */
774 51 : PetscCall(VecDotBegin(Axi,ur,&rAi)); /* r*A*i */
775 51 : PetscCall(VecDotBegin(Axi,ui,&iAi)); /* i*A*i */
776 51 : PetscCall(VecDotBegin(Bxr,ur,&rBr)); /* r*B*r */
777 51 : PetscCall(VecDotBegin(Bxr,ui,&iBr)); /* i*B*r */
778 51 : PetscCall(VecDotBegin(Bxi,ur,&rBi)); /* r*B*i */
779 51 : PetscCall(VecDotBegin(Bxi,ui,&iBi)); /* i*B*i */
780 51 : PetscCall(VecDotEnd(Axr,ur,&rAr)); /* r*A*r */
781 51 : PetscCall(VecDotEnd(Axr,ui,&iAr)); /* i*A*r */
782 51 : PetscCall(VecDotEnd(Axi,ur,&rAi)); /* r*A*i */
783 51 : PetscCall(VecDotEnd(Axi,ui,&iAi)); /* i*A*i */
784 51 : PetscCall(VecDotEnd(Bxr,ur,&rBr)); /* r*B*r */
785 51 : PetscCall(VecDotEnd(Bxr,ui,&iBr)); /* i*B*r */
786 51 : PetscCall(VecDotEnd(Bxi,ur,&rBi)); /* r*B*i */
787 51 : PetscCall(VecDotEnd(Bxi,ui,&iBi)); /* i*B*i */
788 51 : b0 = rAr+iAi; /* rAr+iAi */
789 51 : b2 = rAi-iAr; /* rAi-iAr */
790 51 : b4 = rBr+iBi; /* rBr+iBi */
791 51 : b6 = rBi-iBr; /* rBi-iBr */
792 51 : b7 = b4*b4 + b6*b6; /* k */
793 51 : *eigr = (b0*b4 + b2*b6) / b7; /* eig_r */
794 51 : *eigi = (b2*b4 - b0*b6) / b7; /* eig_i */
795 51 : PetscFunctionReturn(PETSC_SUCCESS);
796 : }
797 : #endif
798 :
799 4826 : static inline PetscErrorCode dvd_compute_n_rr(PetscInt i_s,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,Vec *u,Vec *Ax,Vec *Bx)
800 : {
801 4826 : PetscInt i;
802 4826 : PetscScalar b0,b1;
803 :
804 4826 : PetscFunctionBegin;
805 9652 : for (i=0; i<n; i++) {
806 : #if !defined(PETSC_USE_COMPLEX)
807 4826 : if (eigi[i_s+i] != 0.0) {
808 51 : PetscScalar eigr0=0.0,eigi0=0.0;
809 51 : PetscCall(dvd_complex_rayleigh_quotient(u[i],u[i+1],Ax[i],Ax[i+1],Bx[i],Bx[i+1],&eigr0,&eigi0));
810 51 : if (PetscAbsScalar(eigr[i_s+i]-eigr0)/PetscAbsScalar(eigr[i_s+i]) > 1e-10 || PetscAbsScalar(eigi[i_s+i]-eigi0)/PetscAbsScalar(eigi[i_s+i]) > 1e-10) PetscCall(PetscInfo(u[0],"The eigenvalue %g%+gi is far from its Rayleigh quotient value %g%+gi\n",(double)eigr[i_s+i],(double)eigi[i_s+i],(double)eigr0,(double)eigi0));
811 51 : i++;
812 : } else
813 : #endif
814 : {
815 4775 : PetscCall(VecDotBegin(Ax[i],u[i],&b0));
816 4775 : PetscCall(VecDotBegin(Bx[i],u[i],&b1));
817 4775 : PetscCall(VecDotEnd(Ax[i],u[i],&b0));
818 4775 : PetscCall(VecDotEnd(Bx[i],u[i],&b1));
819 4775 : b0 = b0/b1;
820 4826 : if (PetscAbsScalar(eigr[i_s+i]-b0)/PetscAbsScalar(eigr[i_s+i]) > 1e-10) PetscCall(PetscInfo(u[0],"The eigenvalue %g+%g is far from its Rayleigh quotient value %g+%g\n",(double)PetscRealPart(eigr[i_s+i]),(double)PetscImaginaryPart(eigr[i_s+i]),(double)PetscRealPart(b0),(double)PetscImaginaryPart(b0)));
821 : }
822 : }
823 4826 : PetscFunctionReturn(PETSC_SUCCESS);
824 : }
825 :
826 : /*
827 : Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
828 : kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
829 : where
830 : pX,pY, the right and left eigenvectors of the projected system
831 : ld, the leading dimension of pX and pY
832 : */
833 4826 : static PetscErrorCode dvd_improvex_jd_proj_uv_KZX(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u,Vec *v,Vec *kr,PetscScalar *theta,PetscScalar *thetai,PetscScalar *pX,PetscScalar *pY,PetscInt ld)
834 : {
835 4826 : PetscInt n = i_e-i_s,i;
836 4826 : PetscScalar *b;
837 4826 : Vec *Ax,*Bx,*r;
838 4826 : Mat M;
839 4826 : BV X;
840 :
841 4826 : PetscFunctionBegin;
842 4826 : PetscCall(BVDuplicateResize(d->eps->V,4,&X));
843 4826 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,4,4,NULL,&M));
844 : /* u <- X(i) */
845 4826 : PetscCall(dvd_improvex_compute_X(d,i_s,i_e,u,pX,ld));
846 :
847 : /* v <- theta[0]A*u + theta[1]*B*u */
848 :
849 : /* Bx <- B*X(i) */
850 4826 : Bx = kr;
851 4826 : if (d->BX) {
852 2322 : for (i=i_s; i<i_e; ++i) PetscCall(BVMultVec(d->BX,1.0,0.0,Bx[i-i_s],&pX[ld*i]));
853 : } else {
854 7381 : for (i=0;i<n;i++) {
855 3693 : if (d->B) PetscCall(MatMult(d->B, u[i], Bx[i]));
856 3693 : else PetscCall(VecCopy(u[i], Bx[i]));
857 : }
858 : }
859 :
860 : /* Ax <- A*X(i) */
861 4826 : PetscCall(SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&r));
862 4826 : Ax = r;
863 9703 : for (i=i_s; i<i_e; ++i) PetscCall(BVMultVec(d->AX,1.0,0.0,Ax[i-i_s],&pX[ld*i]));
864 :
865 : /* v <- Y(i) */
866 9703 : for (i=i_s; i<i_e; ++i) PetscCall(BVMultVec(d->W?d->W:d->eps->V,1.0,0.0,v[i-i_s],&pY[ld*i]));
867 :
868 : /* Recompute the eigenvalue */
869 4826 : PetscCall(dvd_compute_n_rr(i_s,n,d->eigr,d->eigi,v,Ax,Bx));
870 :
871 9652 : for (i=0;i<n;i++) {
872 : #if !defined(PETSC_USE_COMPLEX)
873 4826 : if (d->eigi[i_s+i] != 0.0) {
874 : /* [r_i r_i+1 kr_i kr_i+1]*= [ theta_2i' 0 1 0
875 : 0 theta_2i' 0 1
876 : theta_2i+1 -thetai_i -eigr_i -eigi_i
877 : thetai_i theta_2i+1 eigi_i -eigr_i ] */
878 51 : PetscCall(MatDenseGetArrayWrite(M,&b));
879 51 : b[0] = b[5] = PetscConj(theta[2*i]);
880 51 : b[2] = b[7] = -theta[2*i+1];
881 51 : b[6] = -(b[3] = thetai[i]);
882 51 : b[1] = b[4] = 0.0;
883 51 : b[8] = b[13] = 1.0/d->nX[i_s+i];
884 51 : b[10] = b[15] = -d->eigr[i_s+i]/d->nX[i_s+i];
885 51 : b[14] = -(b[11] = d->eigi[i_s+i]/d->nX[i_s+i]);
886 51 : b[9] = b[12] = 0.0;
887 51 : PetscCall(MatDenseRestoreArrayWrite(M,&b));
888 51 : PetscCall(BVInsertVec(X,0,Ax[i]));
889 51 : PetscCall(BVInsertVec(X,1,Ax[i+1]));
890 51 : PetscCall(BVInsertVec(X,2,Bx[i]));
891 51 : PetscCall(BVInsertVec(X,3,Bx[i+1]));
892 51 : PetscCall(BVSetActiveColumns(X,0,4));
893 51 : PetscCall(BVMultInPlace(X,M,0,4));
894 51 : PetscCall(BVCopyVec(X,0,Ax[i]));
895 51 : PetscCall(BVCopyVec(X,1,Ax[i+1]));
896 51 : PetscCall(BVCopyVec(X,2,Bx[i]));
897 51 : PetscCall(BVCopyVec(X,3,Bx[i+1]));
898 51 : i++;
899 : } else
900 : #endif
901 : {
902 : /* [Ax_i Bx_i]*= [ theta_2i' 1/nX_i
903 : theta_2i+1 -eig_i/nX_i ] */
904 4775 : PetscCall(MatDenseGetArrayWrite(M,&b));
905 4775 : b[0] = PetscConj(theta[i*2]);
906 4775 : b[1] = theta[i*2+1];
907 4775 : b[4] = 1.0/d->nX[i_s+i];
908 4775 : b[5] = -d->eigr[i_s+i]/d->nX[i_s+i];
909 4775 : PetscCall(MatDenseRestoreArrayWrite(M,&b));
910 4775 : PetscCall(BVInsertVec(X,0,Ax[i]));
911 4775 : PetscCall(BVInsertVec(X,1,Bx[i]));
912 4775 : PetscCall(BVSetActiveColumns(X,0,2));
913 4775 : PetscCall(BVMultInPlace(X,M,0,2));
914 4775 : PetscCall(BVCopyVec(X,0,Ax[i]));
915 4826 : PetscCall(BVCopyVec(X,1,Bx[i]));
916 : }
917 : }
918 9703 : for (i=0; i<n; i++) d->nX[i_s+i] = 1.0;
919 :
920 : /* v <- K^{-1} r = K^{-1}(theta_2i'*Ax + theta_2i+1*Bx) */
921 9703 : for (i=0;i<n;i++) PetscCall(d->improvex_precond(d,i_s+i,r[i],v[i]));
922 4826 : PetscCall(SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&r));
923 :
924 : /* kr <- P*(Ax - eig_i*Bx) */
925 4826 : PetscCall(d->calcpairs_proj_res(d,i_s,i_e,kr));
926 4826 : PetscCall(BVDestroy(&X));
927 4826 : PetscCall(MatDestroy(&M));
928 4826 : PetscFunctionReturn(PETSC_SUCCESS);
929 : }
930 :
931 4877 : static PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d,PetscInt i,PetscScalar* theta,PetscScalar* thetai,PetscInt *maxits,PetscReal *tol)
932 : {
933 4877 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
934 4877 : PetscReal a;
935 :
936 4877 : PetscFunctionBegin;
937 4877 : a = SlepcAbsEigenvalue(d->eigr[i],d->eigi[i]);
938 :
939 4877 : if (d->nR[i] < data->fix*a) {
940 353 : theta[0] = d->eigr[i];
941 353 : theta[1] = 1.0;
942 : #if !defined(PETSC_USE_COMPLEX)
943 353 : *thetai = d->eigi[i];
944 : #endif
945 : } else {
946 4524 : theta[0] = d->target[0];
947 4524 : theta[1] = d->target[1];
948 : #if !defined(PETSC_USE_COMPLEX)
949 4524 : *thetai = 0.0;
950 : #endif
951 : }
952 :
953 : #if defined(PETSC_USE_COMPLEX)
954 : if (thetai) *thetai = 0.0;
955 : #endif
956 4877 : *maxits = data->maxits;
957 4877 : *tol = data->tol;
958 4877 : PetscFunctionReturn(PETSC_SUCCESS);
959 : }
960 :
961 234 : PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard *d,dvdBlackboard *b,PetscInt maxits,PetscReal tol,PetscReal fix)
962 : {
963 234 : dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
964 :
965 234 : PetscFunctionBegin;
966 : /* Setup the step */
967 234 : if (b->state >= DVD_STATE_CONF) {
968 78 : data->maxits = maxits;
969 78 : data->tol = tol;
970 78 : data->fix = fix;
971 78 : d->improvex_jd_lit = dvd_improvex_jd_lit_const_0;
972 : }
973 234 : PetscFunctionReturn(PETSC_SUCCESS);
974 : }
975 :
976 234 : PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d,dvdBlackboard *b)
977 : {
978 234 : PetscFunctionBegin;
979 : /* Setup the step */
980 234 : if (b->state >= DVD_STATE_CONF) {
981 78 : d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX;
982 : }
983 234 : PetscFunctionReturn(PETSC_SUCCESS);
984 : }
985 :
986 23477 : PetscErrorCode dvd_improvex_compute_X(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u_,PetscScalar *pX,PetscInt ld)
987 : {
988 23477 : PetscInt n = i_e - i_s,i;
989 23477 : Vec *u;
990 :
991 23477 : PetscFunctionBegin;
992 23477 : if (u_) u = u_;
993 2453 : else if (d->correctXnorm) PetscCall(SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&u));
994 23477 : if (u_ || d->correctXnorm) {
995 42383 : for (i=0;i<n;i++) PetscCall(BVMultVec(d->eps->V,1.0,0.0,u[i],&pX[ld*(i+i_s)]));
996 : }
997 : /* nX(i) <- ||X(i)|| */
998 23477 : if (d->correctXnorm) {
999 2434 : for (i=0;i<n;i++) PetscCall(VecNormBegin(u[i],NORM_2,&d->nX[i_s+i]));
1000 2434 : for (i=0;i<n;i++) PetscCall(VecNormEnd(u[i],NORM_2,&d->nX[i_s+i]));
1001 : #if !defined(PETSC_USE_COMPLEX)
1002 2434 : for (i=0;i<n;i++) {
1003 1217 : if (d->eigi[i_s+i] != 0.0) {
1004 0 : d->nX[i_s+i] = d->nX[i_s+i+1] = PetscSqrtScalar(d->nX[i_s+i]*d->nX[i_s+i]+d->nX[i_s+i+1]*d->nX[i_s+i+1]);
1005 0 : i++;
1006 : }
1007 : }
1008 : #endif
1009 : } else {
1010 44591 : for (i=0;i<n;i++) d->nX[i_s+i] = 1.0;
1011 : }
1012 23477 : if (d->correctXnorm && !u_) PetscCall(SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&u));
1013 23477 : PetscFunctionReturn(PETSC_SUCCESS);
1014 : }
|
