Actual source code: dvdimprovex.c
slepc-3.21.1 2024-04-26
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,info,ld;
58: #if defined(PETSC_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 defined(PETSC_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: PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
85: PetscCall(PetscFPTrapPop());
86: SlepcCheckLapackInfo("getrs",info);
88: /* V <- V - KZ*h */
89: PetscCall(BVSetActiveColumns(data->KZ,0,k));
90: for (i=0;i<cV;i++) PetscCall(BVMultVec(data->KZ,-1.0,1.0,V[i],&h[ldh*i]));
91: PetscCall(BVSetActiveColumns(data->KZ,l,k));
92: PetscCall(PetscFree(h));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
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: static PetscErrorCode dvd_improvex_applytrans_proj(dvdDashboard *d,Vec *V,PetscInt cV)
102: {
103: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
104: PetscInt i,ldh,k,l;
105: PetscScalar *h;
106: PetscBLASInt cV_, n, info, ld;
107: #if defined(PETSC_USE_COMPLEX)
108: PetscInt j;
109: #endif
111: PetscFunctionBegin;
112: PetscAssert(cV<=2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
114: /* h <- KZ'*V */
115: PetscCall(PetscMalloc1(data->size_iXKZ*cV,&h));
116: ldh = data->size_iXKZ;
117: PetscCall(BVGetActiveColumns(data->U,&l,&k));
118: PetscAssert(ldh==k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
119: PetscCall(BVSetActiveColumns(data->KZ,0,k));
120: for (i=0;i<cV;i++) {
121: 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: PetscCall(BVSetActiveColumns(data->KZ,l,k));
128: /* h <- iXKZ\h */
129: PetscCall(PetscBLASIntCast(cV,&cV_));
130: PetscCall(PetscBLASIntCast(data->size_iXKZ,&n));
131: PetscCall(PetscBLASIntCast(data->ldiXKZ,&ld));
132: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
133: PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
134: PetscCall(PetscFPTrapPop());
135: SlepcCheckLapackInfo("getrs",info);
137: /* V <- V - U*h */
138: PetscCall(BVSetActiveColumns(data->U,0,k));
139: for (i=0;i<cV;i++) PetscCall(BVMultVec(data->U,-1.0,1.0,V[i],&h[ldh*i]));
140: PetscCall(BVSetActiveColumns(data->U,l,k));
141: PetscCall(PetscFree(h));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)
146: {
147: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
149: PetscFunctionBegin;
150: PetscCall(VecDestroy(&data->friends));
152: /* Restore the pc of ksp */
153: if (data->old_pc) {
154: PetscCall(KSPSetPC(data->ksp, data->old_pc));
155: PetscCall(PCDestroy(&data->old_pc));
156: }
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)
161: {
162: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
164: PetscFunctionBegin;
165: /* Free local data and objects */
166: PetscCall(PetscFree(data->XKZ));
167: PetscCall(PetscFree(data->iXKZ));
168: PetscCall(PetscFree(data->iXKZPivots));
169: PetscCall(BVDestroy(&data->KZ));
170: PetscCall(BVDestroy(&data->U));
171: PetscCall(PetscFree(data));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*
176: y <- theta[1]A*x - theta[0]*B*x
177: auxV, two auxiliary vectors
178: */
179: static inline PetscErrorCode dvd_aux_matmult(dvdImprovex_jd *data,const Vec *x,const Vec *y)
180: {
181: PetscInt n,i;
182: const Vec *Bx;
183: Vec *auxV;
185: PetscFunctionBegin;
186: n = data->r_e - data->r_s;
187: for (i=0;i<n;i++) PetscCall(MatMult(data->d->A,x[i],y[i]));
189: PetscCall(SlepcVecPoolGetVecs(data->d->auxV,2,&auxV));
190: for (i=0;i<n;i++) {
191: #if !defined(PETSC_USE_COMPLEX)
192: if (PetscUnlikely(data->d->eigi[data->r_s+i] != 0.0)) {
193: if (data->d->B) {
194: PetscCall(MatMult(data->d->B,x[i],auxV[0]));
195: PetscCall(MatMult(data->d->B,x[i+1],auxV[1]));
196: Bx = auxV;
197: } else Bx = &x[i];
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: PetscCall(VecAXPBYPCZ(y[i],-data->theta[2*i],data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]));
202: PetscCall(VecAXPBYPCZ(y[i+1],-data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]));
203: i++;
204: } else
205: #endif
206: {
207: if (data->d->B) {
208: PetscCall(MatMult(data->d->B,x[i],auxV[0]));
209: Bx = auxV;
210: } else Bx = &x[i];
211: PetscCall(VecAXPBY(y[i],-data->theta[i*2],data->theta[i*2+1],Bx[0]));
212: }
213: }
214: PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*
219: y <- theta[1]'*A'*x - theta[0]'*B'*x
220: */
221: static inline PetscErrorCode dvd_aux_matmulttrans(dvdImprovex_jd *data,const Vec *x,const Vec *y)
222: {
223: PetscInt n,i;
224: const Vec *Bx;
225: Vec *auxV;
227: PetscFunctionBegin;
228: n = data->r_e - data->r_s;
229: for (i=0;i<n;i++) PetscCall(MatMultTranspose(data->d->A,x[i],y[i]));
231: PetscCall(SlepcVecPoolGetVecs(data->d->auxV,2,&auxV));
232: for (i=0;i<n;i++) {
233: #if !defined(PETSC_USE_COMPLEX)
234: if (data->d->eigi[data->r_s+i] != 0.0) {
235: if (data->d->B) {
236: PetscCall(MatMultTranspose(data->d->B,x[i],auxV[0]));
237: PetscCall(MatMultTranspose(data->d->B,x[i+1],auxV[1]));
238: Bx = auxV;
239: } else Bx = &x[i];
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: PetscCall(VecAXPBYPCZ(y[i],-data->theta[2*i],-data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]));
244: PetscCall(VecAXPBYPCZ(y[i+1],data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]));
245: i++;
246: } else
247: #endif
248: {
249: if (data->d->B) {
250: PetscCall(MatMultTranspose(data->d->B,x[i],auxV[0]));
251: Bx = auxV;
252: } else Bx = &x[i];
253: PetscCall(VecAXPBY(y[i],PetscConj(-data->theta[i*2]),PetscConj(data->theta[i*2+1]),Bx[0]));
254: }
255: }
256: PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode PCApplyBA_dvd(PC pc,PCSide side,Vec in,Vec out,Vec w)
261: {
262: dvdImprovex_jd *data;
263: PetscInt n,i;
264: const Vec *inx,*outx,*wx;
265: Vec *auxV;
266: Mat A;
268: PetscFunctionBegin;
269: PetscCall(PCGetOperators(pc,&A,NULL));
270: PetscCall(MatShellGetContext(A,&data));
271: PetscCall(VecCompGetSubVecs(in,NULL,&inx));
272: PetscCall(VecCompGetSubVecs(out,NULL,&outx));
273: PetscCall(VecCompGetSubVecs(w,NULL,&wx));
274: n = data->r_e - data->r_s;
275: PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
276: switch (side) {
277: case PC_LEFT:
278: /* aux <- theta[1]A*in - theta[0]*B*in */
279: PetscCall(dvd_aux_matmult(data,inx,auxV));
281: /* out <- K * aux */
282: 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: for (i=0;i<n;i++) PetscCall(data->d->improvex_precond(data->d,data->r_s+i,inx[i],auxV[i]));
288: /* out <- theta[1]A*auxV - theta[0]*B*auxV */
289: PetscCall(dvd_aux_matmult(data,auxV,outx));
290: break;
291: case PC_SYMMETRIC:
292: /* aux <- K^{1/2} * in */
293: for (i=0;i<n;i++) PetscCall(PCApplySymmetricRight(data->old_pc,inx[i],auxV[i]));
295: /* wx <- theta[1]A*auxV - theta[0]*B*auxV */
296: PetscCall(dvd_aux_matmult(data,auxV,wx));
298: /* aux <- K^{1/2} * in */
299: for (i=0;i<n;i++) PetscCall(PCApplySymmetricLeft(data->old_pc,wx[i],outx[i]));
300: break;
301: default:
302: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported KSP side");
303: }
304: /* out <- out - v*(u'*out) */
305: PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
306: PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode PCApply_dvd(PC pc,Vec in,Vec out)
311: {
312: dvdImprovex_jd *data;
313: PetscInt n,i;
314: const Vec *inx, *outx;
315: Mat A;
317: PetscFunctionBegin;
318: PetscCall(PCGetOperators(pc,&A,NULL));
319: PetscCall(MatShellGetContext(A,&data));
320: PetscCall(VecCompGetSubVecs(in,NULL,&inx));
321: PetscCall(VecCompGetSubVecs(out,NULL,&outx));
322: n = data->r_e - data->r_s;
323: /* out <- K * in */
324: 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: PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode PCApplyTranspose_dvd(PC pc,Vec in,Vec out)
331: {
332: dvdImprovex_jd *data;
333: PetscInt n,i;
334: const Vec *inx, *outx;
335: Vec *auxV;
336: Mat A;
338: PetscFunctionBegin;
339: PetscCall(PCGetOperators(pc,&A,NULL));
340: PetscCall(MatShellGetContext(A,&data));
341: PetscCall(VecCompGetSubVecs(in,NULL,&inx));
342: PetscCall(VecCompGetSubVecs(out,NULL,&outx));
343: n = data->r_e - data->r_s;
344: PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
345: /* auxV <- in */
346: for (i=0;i<n;i++) PetscCall(VecCopy(inx[i],auxV[i]));
347: /* auxV <- auxV - u*(v'*auxV) */
348: PetscCall(dvd_improvex_applytrans_proj(data->d,auxV,n));
349: /* out <- K' * aux */
350: for (i=0;i<n;i++) PetscCall(PCApplyTranspose(data->old_pc,auxV[i],outx[i]));
351: PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode MatMult_dvd_jd(Mat A,Vec in,Vec out)
356: {
357: dvdImprovex_jd *data;
358: PetscInt n;
359: const Vec *inx, *outx;
360: PCSide side;
362: PetscFunctionBegin;
363: PetscCall(MatShellGetContext(A,&data));
364: PetscCall(VecCompGetSubVecs(in,NULL,&inx));
365: PetscCall(VecCompGetSubVecs(out,NULL,&outx));
366: n = data->r_e - data->r_s;
367: /* out <- theta[1]A*in - theta[0]*B*in */
368: PetscCall(dvd_aux_matmult(data,inx,outx));
369: PetscCall(KSPGetPCSide(data->ksp,&side));
370: if (side == PC_RIGHT) {
371: /* out <- out - v*(u'*out) */
372: PetscCall(dvd_improvex_apply_proj(data->d,(Vec*)outx,n));
373: }
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: static PetscErrorCode MatMultTranspose_dvd_jd(Mat A,Vec in,Vec out)
378: {
379: dvdImprovex_jd *data;
380: PetscInt n,i;
381: const Vec *inx,*outx,*r;
382: Vec *auxV;
383: PCSide side;
385: PetscFunctionBegin;
386: PetscCall(MatShellGetContext(A,&data));
387: PetscCall(VecCompGetSubVecs(in,NULL,&inx));
388: PetscCall(VecCompGetSubVecs(out,NULL,&outx));
389: n = data->r_e - data->r_s;
390: PetscCall(KSPGetPCSide(data->ksp,&side));
391: if (side == PC_RIGHT) {
392: /* auxV <- in */
393: PetscCall(SlepcVecPoolGetVecs(data->d->auxV,n,&auxV));
394: for (i=0;i<n;i++) PetscCall(VecCopy(inx[i],auxV[i]));
395: /* auxV <- auxV - v*(u'*auxV) */
396: PetscCall(dvd_improvex_applytrans_proj(data->d,auxV,n));
397: r = auxV;
398: } else r = inx;
399: /* out <- theta[1]A*r - theta[0]*B*r */
400: PetscCall(dvd_aux_matmulttrans(data,r,outx));
401: if (side == PC_RIGHT) PetscCall(SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode MatCreateVecs_dvd_jd(Mat A,Vec *right,Vec *left)
406: {
407: Vec *r,*l;
408: dvdImprovex_jd *data;
409: PetscInt n,i;
411: PetscFunctionBegin;
412: PetscCall(MatShellGetContext(A,&data));
413: n = data->ksp_max_size;
414: if (right) PetscCall(PetscMalloc1(n,&r));
415: if (left) PetscCall(PetscMalloc1(n,&l));
416: for (i=0;i<n;i++) PetscCall(MatCreateVecs(data->d->A,right?&r[i]:NULL,left?&l[i]:NULL));
417: if (right) {
418: PetscCall(VecCreateCompWithVecs(r,n,data->friends,right));
419: for (i=0;i<n;i++) PetscCall(VecDestroy(&r[i]));
420: }
421: if (left) {
422: PetscCall(VecCreateCompWithVecs(l,n,data->friends,left));
423: for (i=0;i<n;i++) PetscCall(VecDestroy(&l[i]));
424: }
426: if (right) PetscCall(PetscFree(r));
427: if (left) PetscCall(PetscFree(l));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)
432: {
433: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
434: PetscInt rA, cA, rlA, clA;
435: Mat A;
436: PetscBool t;
437: PC pc;
438: Vec v0[2];
440: PetscFunctionBegin;
441: data->size_cX = data->old_size_X = 0;
442: data->lastTol = data->dynamic?0.5:0.0;
444: /* Setup the ksp */
445: if (data->ksp) {
446: /* Create the reference vector */
447: PetscCall(BVGetColumn(d->eps->V,0,&v0[0]));
448: v0[1] = v0[0];
449: PetscCall(VecCreateCompWithVecs(v0,data->ksp_max_size,NULL,&data->friends));
450: PetscCall(BVRestoreColumn(d->eps->V,0,&v0[0]));
452: /* Save the current pc and set a PCNONE */
453: PetscCall(KSPGetPC(data->ksp, &data->old_pc));
454: PetscCall(PetscObjectTypeCompare((PetscObject)data->old_pc,PCNONE,&t));
455: data->lastTol = 0.5;
456: if (t) data->old_pc = NULL;
457: else {
458: PetscCall(PetscObjectReference((PetscObject)data->old_pc));
459: PetscCall(PCCreate(PetscObjectComm((PetscObject)d->eps),&pc));
460: PetscCall(PCSetType(pc,PCSHELL));
461: PetscCall(PCSetOperators(pc,d->A,d->A));
462: PetscCall(PCSetReusePreconditioner(pc,PETSC_TRUE));
463: PetscCall(PCShellSetApply(pc,PCApply_dvd));
464: PetscCall(PCShellSetApplyBA(pc,PCApplyBA_dvd));
465: PetscCall(PCShellSetApplyTranspose(pc,PCApplyTranspose_dvd));
466: PetscCall(KSPSetPC(data->ksp,pc));
467: PetscCall(PCDestroy(&pc));
468: }
470: /* Create the (I-v*u')*K*(A-s*B) matrix */
471: PetscCall(MatGetSize(d->A,&rA,&cA));
472: PetscCall(MatGetLocalSize(d->A,&rlA,&clA));
473: 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: PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_dvd_jd));
475: PetscCall(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_dvd_jd));
476: PetscCall(MatShellSetOperation(A,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_dvd_jd));
478: /* Try to avoid KSPReset */
479: PetscCall(KSPGetOperatorsSet(data->ksp,&t,NULL));
480: if (t) {
481: Mat M;
482: PetscInt rM;
483: PetscCall(KSPGetOperators(data->ksp,&M,NULL));
484: PetscCall(MatGetSize(M,&rM,NULL));
485: if (rM != rA*data->ksp_max_size) PetscCall(KSPReset(data->ksp));
486: }
487: PetscCall(EPS_KSPSetOperators(data->ksp,A,A));
488: PetscCall(KSPSetReusePreconditioner(data->ksp,PETSC_TRUE));
489: PetscCall(KSPSetUp(data->ksp));
490: PetscCall(MatDestroy(&A));
491: } else {
492: data->old_pc = NULL;
493: data->friends = NULL;
494: }
495: PetscCall(BVSetActiveColumns(data->KZ,0,0));
496: PetscCall(BVSetActiveColumns(data->U,0,0));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
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: 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: PetscInt n=i_e-i_s,size_KZ,V_new,rm,i,lv,kv,lKZ,kKZ;
510: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
511: const PetscScalar *array;
512: Mat M;
513: Vec u[2],v[2];
514: PetscBLASInt s,ldXKZ,info;
516: PetscFunctionBegin;
517: /* Check consistency */
518: PetscCall(BVGetActiveColumns(d->eps->V,&lv,&kv));
519: V_new = lv - data->size_cX;
520: PetscAssert(V_new<=data->old_size_X,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
521: data->old_size_X = n;
522: data->size_cX = lv;
524: /* KZ <- KZ(rm:rm+max_cX-1) */
525: PetscCall(BVGetActiveColumns(data->KZ,&lKZ,&kKZ));
526: rm = PetscMax(V_new+lKZ,0);
527: if (rm > 0) {
528: for (i=0;i<lKZ;i++) {
529: PetscCall(BVCopyColumn(data->KZ,i+rm,i));
530: PetscCall(BVCopyColumn(data->U,i+rm,i));
531: }
532: }
534: /* XKZ <- XKZ(rm:rm+max_cX-1,rm:rm+max_cX-1) */
535: if (rm > 0) {
536: 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: lKZ = PetscMin(0,lKZ+V_new);
539: PetscCall(BVSetActiveColumns(data->KZ,lKZ,lKZ+n));
540: PetscCall(BVSetActiveColumns(data->U,lKZ,lKZ+n));
542: /* Compute X, KZ and KR */
543: PetscCall(BVGetColumn(data->U,lKZ,u));
544: if (n>1) PetscCall(BVGetColumn(data->U,lKZ+1,&u[1]));
545: PetscCall(BVGetColumn(data->KZ,lKZ,v));
546: if (n>1) PetscCall(BVGetColumn(data->KZ,lKZ+1,&v[1]));
547: PetscCall(d->improvex_jd_proj_uv(d,i_s,i_e,u,v,kr,theta,thetai,pX,pY,ld));
548: PetscCall(BVRestoreColumn(data->U,lKZ,u));
549: if (n>1) PetscCall(BVRestoreColumn(data->U,lKZ+1,&u[1]));
550: PetscCall(BVRestoreColumn(data->KZ,lKZ,v));
551: if (n>1) PetscCall(BVRestoreColumn(data->KZ,lKZ+1,&v[1]));
553: /* XKZ <- U'*KZ */
554: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,lKZ+n,lKZ+n,NULL,&M));
555: PetscCall(BVMatProject(data->KZ,NULL,data->U,M));
556: PetscCall(MatDenseGetArrayRead(M,&array));
557: for (i=lKZ;i<lKZ+n;i++) { /* upper part */
558: PetscCall(PetscArraycpy(&data->XKZ[data->ldXKZ*i],&array[i*(lKZ+n)],lKZ));
559: }
560: for (i=0;i<lKZ+n;i++) { /* lower part */
561: PetscCall(PetscArraycpy(&data->XKZ[data->ldXKZ*i+lKZ],&array[i*(lKZ+n)+lKZ],n));
562: }
563: PetscCall(MatDenseRestoreArrayRead(M,&array));
564: PetscCall(MatDestroy(&M));
566: /* iXKZ <- inv(XKZ) */
567: size_KZ = lKZ+n;
568: PetscCall(PetscBLASIntCast(lKZ+n,&s));
569: data->ldiXKZ = data->size_iXKZ = size_KZ;
570: for (i=0;i<size_KZ;i++) PetscCall(PetscArraycpy(&data->iXKZ[data->ldiXKZ*i],&data->XKZ[data->ldXKZ*i],size_KZ));
571: PetscCall(PetscBLASIntCast(data->ldiXKZ,&ldXKZ));
572: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
573: PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&s,&s,data->iXKZ,&ldXKZ,data->iXKZPivots,&info));
574: PetscCall(PetscFPTrapPop());
575: SlepcCheckLapackInfo("getrf",info);
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
580: {
581: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
582: PetscInt i,j,n,maxits,maxits0,lits,s,ld,k,max_size_D,lV,kV;
583: PetscScalar *pX,*pY;
584: PetscReal tol,tol0;
585: Vec *kr,kr_comp,D_comp,D[2],kr0[2];
586: PetscBool odd_situation = PETSC_FALSE;
588: PetscFunctionBegin;
589: PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
590: max_size_D = d->eps->ncv-kV;
591: /* Quick exit */
592: if ((max_size_D == 0) || r_e-r_s <= 0) {
593: *size_D = 0;
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);
598: PetscAssert(n>0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"n == 0");
599: PetscAssert(data->size_X>=r_e-r_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"size_X < r_e-r_s");
601: PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
603: /* Restart lastTol if a new pair converged */
604: if (data->dynamic && data->size_cX < lV)
605: data->lastTol = 0.5;
607: 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: if (d->eigi[r_s+i] != 0.0) {
611: if (i+2 <= max_size_D) s=2;
612: else break;
613: } else
614: #endif
615: s=1;
617: data->r_s = r_s+i;
618: data->r_e = r_s+i+s;
619: PetscCall(SlepcVecPoolGetVecs(d->auxV,s,&kr));
621: /* Compute theta, maximum iterations and tolerance */
622: maxits = 0;
623: tol = 1;
624: for (j=0;j<s;j++) {
625: PetscCall(d->improvex_jd_lit(d,r_s+i+j,&data->theta[2*j],&data->thetai[j],&maxits0,&tol0));
626: maxits += maxits0;
627: tol *= tol0;
628: }
629: maxits/= s;
630: tol = data->dynamic?data->lastTol:PetscExpReal(PetscLogReal(tol)/s);
632: /* Compute u, v and kr */
633: k = r_s+i;
634: PetscCall(DSVectors(d->eps->ds,DS_MAT_X,&k,NULL));
635: k = r_s+i;
636: PetscCall(DSVectors(d->eps->ds,DS_MAT_Y,&k,NULL));
637: PetscCall(DSGetArray(d->eps->ds,DS_MAT_X,&pX));
638: PetscCall(DSGetArray(d->eps->ds,DS_MAT_Y,&pY));
639: PetscCall(dvd_improvex_jd_proj_cuv(d,r_s+i,r_s+i+s,kr,data->theta,data->thetai,pX,pY,ld));
640: PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_X,&pX));
641: PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_Y,&pY));
643: /* Check if the first eigenpairs are converged */
644: if (i == 0) {
645: PetscInt oldnpreconv = d->npreconv;
646: PetscCall(d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv));
647: if (d->npreconv > oldnpreconv) break;
648: }
650: /* Test the odd situation of solving Ax=b with A=I */
651: #if !defined(PETSC_USE_COMPLEX)
652: 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: if (data->ksp && !odd_situation) {
658: /* kr <- -kr */
659: for (j=0;j<s;j++) PetscCall(VecScale(kr[j],-1.0));
661: /* Compose kr and D */
662: kr0[0] = kr[0];
663: kr0[1] = (s==2 ? kr[1] : NULL);
664: PetscCall(VecCreateCompWithVecs(kr0,data->ksp_max_size,data->friends,&kr_comp));
665: PetscCall(BVGetColumn(d->eps->V,kV+i,&D[0]));
666: if (s==2) PetscCall(BVGetColumn(d->eps->V,kV+i+1,&D[1]));
667: else D[1] = NULL;
668: PetscCall(VecCreateCompWithVecs(D,data->ksp_max_size,data->friends,&D_comp));
669: PetscCall(VecCompSetSubVecs(data->friends,s,NULL));
671: /* Solve the correction equation */
672: PetscCall(KSPSetTolerances(data->ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,maxits));
673: PetscCall(KSPSolve(data->ksp,kr_comp,D_comp));
674: PetscCall(KSPGetIterationNumber(data->ksp,&lits));
676: /* Destroy the composed ks and D */
677: PetscCall(VecDestroy(&kr_comp));
678: PetscCall(VecDestroy(&D_comp));
679: PetscCall(BVRestoreColumn(d->eps->V,kV+i,&D[0]));
680: if (s==2) PetscCall(BVRestoreColumn(d->eps->V,kV+i+1,&D[1]));
682: /* If GD */
683: } else {
684: PetscCall(BVGetColumn(d->eps->V,kV+i,&D[0]));
685: if (s==2) PetscCall(BVGetColumn(d->eps->V,kV+i+1,&D[1]));
686: for (j=0;j<s;j++) PetscCall(d->improvex_precond(d,r_s+i+j,kr[j],D[j]));
687: PetscCall(dvd_improvex_apply_proj(d,D,s));
688: PetscCall(BVRestoreColumn(d->eps->V,kV+i,&D[0]));
689: 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: if (i == 0 && d->eps->errest[d->nconv+r_s] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s] < PETSC_MAX_REAL) {
693: 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: PetscCall(SlepcVecPoolRestoreVecs(d->auxV,s,&kr));
696: }
697: *size_D = i;
698: if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: PetscErrorCode dvd_improvex_jd(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs,PetscBool dynamic)
703: {
704: dvdImprovex_jd *data;
705: PetscBool useGD;
706: PC pc;
707: PetscInt size_P;
709: PetscFunctionBegin;
710: /* Setting configuration constrains */
711: PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&useGD));
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: if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) {
717: max_bs++;
718: b->max_size_P = PetscMax(b->max_size_P,2);
719: } else
720: #endif
721: {
722: b->max_size_P = PetscMax(b->max_size_P,1);
723: }
724: b->max_size_X = PetscMax(b->max_size_X,max_bs);
725: size_P = b->max_size_P;
727: /* Setup the preconditioner */
728: if (ksp) {
729: PetscCall(KSPGetPC(ksp,&pc));
730: PetscCall(dvd_static_precond_PC(d,b,pc));
731: } else PetscCall(dvd_static_precond_PC(d,b,NULL));
733: /* Setup the step */
734: if (b->state >= DVD_STATE_CONF) {
735: PetscCall(PetscNew(&data));
736: data->dynamic = dynamic;
737: PetscCall(PetscMalloc1(size_P*size_P,&data->XKZ));
738: PetscCall(PetscMalloc1(size_P*size_P,&data->iXKZ));
739: PetscCall(PetscMalloc1(size_P,&data->iXKZPivots));
740: data->ldXKZ = size_P;
741: data->size_X = b->max_size_X;
742: d->improveX_data = data;
743: data->ksp = useGD? NULL: ksp;
744: data->d = d;
745: d->improveX = dvd_improvex_jd_gen;
746: #if !defined(PETSC_USE_COMPLEX)
747: if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) data->ksp_max_size = 2;
748: else
749: #endif
750: data->ksp_max_size = 1;
751: /* Create various vector basis */
752: PetscCall(BVDuplicateResize(d->eps->V,size_P,&data->KZ));
753: PetscCall(BVSetMatrix(data->KZ,NULL,PETSC_FALSE));
754: PetscCall(BVDuplicate(data->KZ,&data->U));
756: PetscCall(EPSDavidsonFLAdd(&d->startList,dvd_improvex_jd_start));
757: PetscCall(EPSDavidsonFLAdd(&d->endList,dvd_improvex_jd_end));
758: PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_jd_d));
759: }
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: #if !defined(PETSC_USE_COMPLEX)
764: 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: PetscScalar rAr,iAr,rAi,iAi,rBr,iBr,rBi,iBi,b0,b2,b4,b6,b7;
768: 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: PetscCall(VecDotBegin(Axr,ur,&rAr)); /* r*A*r */
773: PetscCall(VecDotBegin(Axr,ui,&iAr)); /* i*A*r */
774: PetscCall(VecDotBegin(Axi,ur,&rAi)); /* r*A*i */
775: PetscCall(VecDotBegin(Axi,ui,&iAi)); /* i*A*i */
776: PetscCall(VecDotBegin(Bxr,ur,&rBr)); /* r*B*r */
777: PetscCall(VecDotBegin(Bxr,ui,&iBr)); /* i*B*r */
778: PetscCall(VecDotBegin(Bxi,ur,&rBi)); /* r*B*i */
779: PetscCall(VecDotBegin(Bxi,ui,&iBi)); /* i*B*i */
780: PetscCall(VecDotEnd(Axr,ur,&rAr)); /* r*A*r */
781: PetscCall(VecDotEnd(Axr,ui,&iAr)); /* i*A*r */
782: PetscCall(VecDotEnd(Axi,ur,&rAi)); /* r*A*i */
783: PetscCall(VecDotEnd(Axi,ui,&iAi)); /* i*A*i */
784: PetscCall(VecDotEnd(Bxr,ur,&rBr)); /* r*B*r */
785: PetscCall(VecDotEnd(Bxr,ui,&iBr)); /* i*B*r */
786: PetscCall(VecDotEnd(Bxi,ur,&rBi)); /* r*B*i */
787: PetscCall(VecDotEnd(Bxi,ui,&iBi)); /* i*B*i */
788: b0 = rAr+iAi; /* rAr+iAi */
789: b2 = rAi-iAr; /* rAi-iAr */
790: b4 = rBr+iBi; /* rBr+iBi */
791: b6 = rBi-iBr; /* rBi-iBr */
792: b7 = b4*b4 + b6*b6; /* k */
793: *eigr = (b0*b4 + b2*b6) / b7; /* eig_r */
794: *eigi = (b2*b4 - b0*b6) / b7; /* eig_i */
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
797: #endif
799: static inline PetscErrorCode dvd_compute_n_rr(PetscInt i_s,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,Vec *u,Vec *Ax,Vec *Bx)
800: {
801: PetscInt i;
802: PetscScalar b0,b1;
804: PetscFunctionBegin;
805: for (i=0; i<n; i++) {
806: #if !defined(PETSC_USE_COMPLEX)
807: if (eigi[i_s+i] != 0.0) {
808: PetscScalar eigr0=0.0,eigi0=0.0;
809: PetscCall(dvd_complex_rayleigh_quotient(u[i],u[i+1],Ax[i],Ax[i+1],Bx[i],Bx[i+1],&eigr0,&eigi0));
810: 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: i++;
812: } else
813: #endif
814: {
815: PetscCall(VecDotBegin(Ax[i],u[i],&b0));
816: PetscCall(VecDotBegin(Bx[i],u[i],&b1));
817: PetscCall(VecDotEnd(Ax[i],u[i],&b0));
818: PetscCall(VecDotEnd(Bx[i],u[i],&b1));
819: b0 = b0/b1;
820: 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: PetscFunctionReturn(PETSC_SUCCESS);
824: }
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: 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: PetscInt n = i_e-i_s,i;
836: PetscScalar *b;
837: Vec *Ax,*Bx,*r;
838: Mat M;
839: BV X;
841: PetscFunctionBegin;
842: PetscCall(BVDuplicateResize(d->eps->V,4,&X));
843: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,4,4,NULL,&M));
844: /* u <- X(i) */
845: PetscCall(dvd_improvex_compute_X(d,i_s,i_e,u,pX,ld));
847: /* v <- theta[0]A*u + theta[1]*B*u */
849: /* Bx <- B*X(i) */
850: Bx = kr;
851: if (d->BX) {
852: 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: for (i=0;i<n;i++) {
855: if (d->B) PetscCall(MatMult(d->B, u[i], Bx[i]));
856: else PetscCall(VecCopy(u[i], Bx[i]));
857: }
858: }
860: /* Ax <- A*X(i) */
861: PetscCall(SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&r));
862: Ax = r;
863: for (i=i_s; i<i_e; ++i) PetscCall(BVMultVec(d->AX,1.0,0.0,Ax[i-i_s],&pX[ld*i]));
865: /* v <- Y(i) */
866: 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]));
868: /* Recompute the eigenvalue */
869: PetscCall(dvd_compute_n_rr(i_s,n,d->eigr,d->eigi,v,Ax,Bx));
871: for (i=0;i<n;i++) {
872: #if !defined(PETSC_USE_COMPLEX)
873: 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: PetscCall(MatDenseGetArrayWrite(M,&b));
879: b[0] = b[5] = PetscConj(theta[2*i]);
880: b[2] = b[7] = -theta[2*i+1];
881: b[6] = -(b[3] = thetai[i]);
882: b[1] = b[4] = 0.0;
883: b[8] = b[13] = 1.0/d->nX[i_s+i];
884: b[10] = b[15] = -d->eigr[i_s+i]/d->nX[i_s+i];
885: b[14] = -(b[11] = d->eigi[i_s+i]/d->nX[i_s+i]);
886: b[9] = b[12] = 0.0;
887: PetscCall(MatDenseRestoreArrayWrite(M,&b));
888: PetscCall(BVInsertVec(X,0,Ax[i]));
889: PetscCall(BVInsertVec(X,1,Ax[i+1]));
890: PetscCall(BVInsertVec(X,2,Bx[i]));
891: PetscCall(BVInsertVec(X,3,Bx[i+1]));
892: PetscCall(BVSetActiveColumns(X,0,4));
893: PetscCall(BVMultInPlace(X,M,0,4));
894: PetscCall(BVCopyVec(X,0,Ax[i]));
895: PetscCall(BVCopyVec(X,1,Ax[i+1]));
896: PetscCall(BVCopyVec(X,2,Bx[i]));
897: PetscCall(BVCopyVec(X,3,Bx[i+1]));
898: 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: PetscCall(MatDenseGetArrayWrite(M,&b));
905: b[0] = PetscConj(theta[i*2]);
906: b[1] = theta[i*2+1];
907: b[4] = 1.0/d->nX[i_s+i];
908: b[5] = -d->eigr[i_s+i]/d->nX[i_s+i];
909: PetscCall(MatDenseRestoreArrayWrite(M,&b));
910: PetscCall(BVInsertVec(X,0,Ax[i]));
911: PetscCall(BVInsertVec(X,1,Bx[i]));
912: PetscCall(BVSetActiveColumns(X,0,2));
913: PetscCall(BVMultInPlace(X,M,0,2));
914: PetscCall(BVCopyVec(X,0,Ax[i]));
915: PetscCall(BVCopyVec(X,1,Bx[i]));
916: }
917: }
918: for (i=0; i<n; i++) d->nX[i_s+i] = 1.0;
920: /* v <- K^{-1} r = K^{-1}(theta_2i'*Ax + theta_2i+1*Bx) */
921: for (i=0;i<n;i++) PetscCall(d->improvex_precond(d,i_s+i,r[i],v[i]));
922: PetscCall(SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&r));
924: /* kr <- P*(Ax - eig_i*Bx) */
925: PetscCall(d->calcpairs_proj_res(d,i_s,i_e,kr));
926: PetscCall(BVDestroy(&X));
927: PetscCall(MatDestroy(&M));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: static PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d,PetscInt i,PetscScalar* theta,PetscScalar* thetai,PetscInt *maxits,PetscReal *tol)
932: {
933: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
934: PetscReal a;
936: PetscFunctionBegin;
937: a = SlepcAbsEigenvalue(d->eigr[i],d->eigi[i]);
939: if (d->nR[i] < data->fix*a) {
940: theta[0] = d->eigr[i];
941: theta[1] = 1.0;
942: #if !defined(PETSC_USE_COMPLEX)
943: *thetai = d->eigi[i];
944: #endif
945: } else {
946: theta[0] = d->target[0];
947: theta[1] = d->target[1];
948: #if !defined(PETSC_USE_COMPLEX)
949: *thetai = 0.0;
950: #endif
951: }
953: #if defined(PETSC_USE_COMPLEX)
954: if (thetai) *thetai = 0.0;
955: #endif
956: *maxits = data->maxits;
957: *tol = data->tol;
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard *d,dvdBlackboard *b,PetscInt maxits,PetscReal tol,PetscReal fix)
962: {
963: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
965: PetscFunctionBegin;
966: /* Setup the step */
967: if (b->state >= DVD_STATE_CONF) {
968: data->maxits = maxits;
969: data->tol = tol;
970: data->fix = fix;
971: d->improvex_jd_lit = dvd_improvex_jd_lit_const_0;
972: }
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d,dvdBlackboard *b)
977: {
978: PetscFunctionBegin;
979: /* Setup the step */
980: if (b->state >= DVD_STATE_CONF) {
981: d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX;
982: }
983: PetscFunctionReturn(PETSC_SUCCESS);
984: }
986: PetscErrorCode dvd_improvex_compute_X(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u_,PetscScalar *pX,PetscInt ld)
987: {
988: PetscInt n = i_e - i_s,i;
989: Vec *u;
991: PetscFunctionBegin;
992: if (u_) u = u_;
993: else if (d->correctXnorm) PetscCall(SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&u));
994: if (u_ || d->correctXnorm) {
995: 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: if (d->correctXnorm) {
999: for (i=0;i<n;i++) PetscCall(VecNormBegin(u[i],NORM_2,&d->nX[i_s+i]));
1000: for (i=0;i<n;i++) PetscCall(VecNormEnd(u[i],NORM_2,&d->nX[i_s+i]));
1001: #if !defined(PETSC_USE_COMPLEX)
1002: for (i=0;i<n;i++) {
1003: if (d->eigi[i_s+i] != 0.0) {
1004: 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: i++;
1006: }
1007: }
1008: #endif
1009: } else {
1010: for (i=0;i<n;i++) d->nX[i_s+i] = 1.0;
1011: }
1012: if (d->correctXnorm && !u_) PetscCall(SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&u));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }