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