Actual source code: dvdimprovex.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 currect 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: #if defined(PETSC_MISSING_LAPACK_GETRS)
56: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable");
57: #else
59: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
60: PetscInt i,ldh,k,l;
61: PetscScalar *h;
62: PetscBLASInt cV_,n,info,ld;
63: #if defined(PETSC_USE_COMPLEX)
64: PetscInt j;
65: #endif
68: if (cV > 2) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
70: /* h <- X'*V */
71: PetscMalloc1(data->size_iXKZ*cV,&h);
72: ldh = data->size_iXKZ;
73: BVGetActiveColumns(data->U,&l,&k);
74: if (ldh!=k) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
75: BVSetActiveColumns(data->U,0,k);
76: for (i=0;i<cV;i++) {
77: BVDotVec(data->U,V[i],&h[ldh*i]);
78: #if defined(PETSC_USE_COMPLEX)
79: for (j=0; j<k; j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
80: #endif
81: }
82: BVSetActiveColumns(data->U,l,k);
84: /* h <- iXKZ\h */
85: PetscBLASIntCast(cV,&cV_);
86: PetscBLASIntCast(data->size_iXKZ,&n);
87: PetscBLASIntCast(data->ldiXKZ,&ld);
89: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
90: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
91: PetscFPTrapPop();
92: SlepcCheckLapackInfo("getrs",info);
94: /* V <- V - KZ*h */
95: BVSetActiveColumns(data->KZ,0,k);
96: for (i=0;i<cV;i++) {
97: BVMultVec(data->KZ,-1.0,1.0,V[i],&h[ldh*i]);
98: }
99: BVSetActiveColumns(data->KZ,l,k);
100: PetscFree(h);
101: return(0);
102: #endif
103: }
105: /*
106: Compute (I - X*iXKZ*KZ')*V where,
107: V, the vectors to apply the projector,
108: cV, the number of vectors in V,
109: */
110: static PetscErrorCode dvd_improvex_applytrans_proj(dvdDashboard *d,Vec *V,PetscInt cV)
111: {
112: #if defined(PETSC_MISSING_LAPACK_GETRS)
114: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable");
115: return(0);
116: #else
118: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
119: PetscInt i,ldh,k,l;
120: PetscScalar *h;
121: PetscBLASInt cV_, n, info, ld;
122: #if defined(PETSC_USE_COMPLEX)
123: PetscInt j;
124: #endif
127: if (cV > 2) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
129: /* h <- KZ'*V */
130: PetscMalloc1(data->size_iXKZ*cV,&h);
131: ldh = data->size_iXKZ;
132: BVGetActiveColumns(data->U,&l,&k);
133: if (ldh!=k) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
134: BVSetActiveColumns(data->KZ,0,k);
135: for (i=0;i<cV;i++) {
136: BVDotVec(data->KZ,V[i],&h[ldh*i]);
137: #if defined(PETSC_USE_COMPLEX)
138: for (j=0;j<k;j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
139: #endif
140: }
141: BVSetActiveColumns(data->KZ,l,k);
143: /* h <- iXKZ\h */
144: PetscBLASIntCast(cV,&cV_);
145: PetscBLASIntCast(data->size_iXKZ,&n);
146: PetscBLASIntCast(data->ldiXKZ,&ld);
148: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
149: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
150: PetscFPTrapPop();
151: SlepcCheckLapackInfo("getrs",info);
153: /* V <- V - U*h */
154: BVSetActiveColumns(data->U,0,k);
155: for (i=0;i<cV;i++) {
156: BVMultVec(data->U,-1.0,1.0,V[i],&h[ldh*i]);
157: }
158: BVSetActiveColumns(data->U,l,k);
159: PetscFree(h);
160: return(0);
161: #endif
162: }
164: static PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)
165: {
167: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
170: VecDestroy(&data->friends);
172: /* Restore the pc of ksp */
173: if (data->old_pc) {
174: KSPSetPC(data->ksp, data->old_pc);
175: PCDestroy(&data->old_pc);
176: }
177: return(0);
178: }
180: static PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)
181: {
183: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
186: /* Free local data and objects */
187: PetscFree(data->XKZ);
188: PetscFree(data->iXKZ);
189: PetscFree(data->iXKZPivots);
190: BVDestroy(&data->KZ);
191: BVDestroy(&data->U);
192: PetscFree(data);
193: return(0);
194: }
196: /*
197: y <- theta[1]A*x - theta[0]*B*x
198: auxV, two auxiliary vectors
199: */
200: PETSC_STATIC_INLINE PetscErrorCode dvd_aux_matmult(dvdImprovex_jd *data,const Vec *x,const Vec *y)
201: {
203: PetscInt n,i;
204: const Vec *Bx;
205: Vec *auxV;
208: n = data->r_e - data->r_s;
209: for (i=0;i<n;i++) {
210: MatMult(data->d->A,x[i],y[i]);
211: }
213: SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);
214: for (i=0;i<n;i++) {
215: #if !defined(PETSC_USE_COMPLEX)
216: if (data->d->eigi[data->r_s+i] != 0.0) {
217: if (data->d->B) {
218: MatMult(data->d->B,x[i],auxV[0]);
219: MatMult(data->d->B,x[i+1],auxV[1]);
220: Bx = auxV;
221: } else Bx = &x[i];
223: /* y_i <- [ t_2i+1*A*x_i - t_2i*Bx_i + ti_i*Bx_i+1;
224: y_i+1 t_2i+1*A*x_i+1 - ti_i*Bx_i - t_2i*Bx_i+1 ] */
225: VecAXPBYPCZ(y[i],-data->theta[2*i],data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);
226: VecAXPBYPCZ(y[i+1],-data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);
227: i++;
228: } else
229: #endif
230: {
231: if (data->d->B) {
232: MatMult(data->d->B,x[i],auxV[0]);
233: Bx = auxV;
234: } else Bx = &x[i];
235: VecAXPBY(y[i],-data->theta[i*2],data->theta[i*2+1],Bx[0]);
236: }
237: }
238: SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);
239: return(0);
240: }
242: /*
243: y <- theta[1]'*A'*x - theta[0]'*B'*x
244: */
245: PETSC_STATIC_INLINE PetscErrorCode dvd_aux_matmulttrans(dvdImprovex_jd *data,const Vec *x,const Vec *y)
246: {
248: PetscInt n,i;
249: const Vec *Bx;
250: Vec *auxV;
253: n = data->r_e - data->r_s;
254: for (i=0;i<n;i++) {
255: MatMultTranspose(data->d->A,x[i],y[i]);
256: }
258: SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);
259: for (i=0;i<n;i++) {
260: #if !defined(PETSC_USE_COMPLEX)
261: if (data->d->eigi[data->r_s+i] != 0.0) {
262: if (data->d->B) {
263: MatMultTranspose(data->d->B,x[i],auxV[0]);
264: MatMultTranspose(data->d->B,x[i+1],auxV[1]);
265: Bx = auxV;
266: } else Bx = &x[i];
268: /* y_i <- [ t_2i+1*A*x_i - t_2i*Bx_i - ti_i*Bx_i+1;
269: y_i+1 t_2i+1*A*x_i+1 + ti_i*Bx_i - t_2i*Bx_i+1 ] */
270: VecAXPBYPCZ(y[i],-data->theta[2*i],-data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);
271: VecAXPBYPCZ(y[i+1],data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);
272: i++;
273: } else
274: #endif
275: {
276: if (data->d->B) {
277: MatMultTranspose(data->d->B,x[i],auxV[0]);
278: Bx = auxV;
279: } else Bx = &x[i];
280: VecAXPBY(y[i],PetscConj(-data->theta[i*2]),PetscConj(data->theta[i*2+1]),Bx[0]);
281: }
282: }
283: SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);
284: return(0);
285: }
287: static PetscErrorCode PCApplyBA_dvd(PC pc,PCSide side,Vec in,Vec out,Vec w)
288: {
290: dvdImprovex_jd *data;
291: PetscInt n,i;
292: const Vec *inx,*outx,*wx;
293: Vec *auxV;
294: Mat A;
297: PCGetOperators(pc,&A,NULL);
298: MatShellGetContext(A,(void**)&data);
299: VecCompGetSubVecs(in,NULL,&inx);
300: VecCompGetSubVecs(out,NULL,&outx);
301: VecCompGetSubVecs(w,NULL,&wx);
302: n = data->r_e - data->r_s;
303: SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
304: switch (side) {
305: case PC_LEFT:
306: /* aux <- theta[1]A*in - theta[0]*B*in */
307: dvd_aux_matmult(data,inx,auxV);
309: /* out <- K * aux */
310: for (i=0;i<n;i++) {
311: data->d->improvex_precond(data->d,data->r_s+i,auxV[i],outx[i]);
312: }
313: break;
314: case PC_RIGHT:
315: /* aux <- K * in */
316: for (i=0;i<n;i++) {
317: data->d->improvex_precond(data->d,data->r_s+i,inx[i],auxV[i]);
318: }
320: /* out <- theta[1]A*auxV - theta[0]*B*auxV */
321: dvd_aux_matmult(data,auxV,outx);
322: break;
323: case PC_SYMMETRIC:
324: /* aux <- K^{1/2} * in */
325: for (i=0;i<n;i++) {
326: PCApplySymmetricRight(data->old_pc,inx[i],auxV[i]);
327: }
329: /* wx <- theta[1]A*auxV - theta[0]*B*auxV */
330: dvd_aux_matmult(data,auxV,wx);
332: /* aux <- K^{1/2} * in */
333: for (i=0;i<n;i++) {
334: PCApplySymmetricLeft(data->old_pc,wx[i],outx[i]);
335: }
336: break;
337: default:
338: SETERRQ(PETSC_COMM_SELF,1,"Unsupported KSP side");
339: }
340: /* out <- out - v*(u'*out) */
341: dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
342: SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
343: return(0);
344: }
346: static PetscErrorCode PCApply_dvd(PC pc,Vec in,Vec out)
347: {
349: dvdImprovex_jd *data;
350: PetscInt n,i;
351: const Vec *inx, *outx;
352: Mat A;
355: PCGetOperators(pc,&A,NULL);
356: MatShellGetContext(A,(void**)&data);
357: VecCompGetSubVecs(in,NULL,&inx);
358: VecCompGetSubVecs(out,NULL,&outx);
359: n = data->r_e - data->r_s;
360: /* out <- K * in */
361: for (i=0;i<n;i++) {
362: data->d->improvex_precond(data->d,data->r_s+i,inx[i],outx[i]);
363: }
364: /* out <- out - v*(u'*out) */
365: dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
366: return(0);
367: }
369: static PetscErrorCode PCApplyTranspose_dvd(PC pc,Vec in,Vec out)
370: {
372: dvdImprovex_jd *data;
373: PetscInt n,i;
374: const Vec *inx, *outx;
375: Vec *auxV;
376: Mat A;
379: PCGetOperators(pc,&A,NULL);
380: MatShellGetContext(A,(void**)&data);
381: VecCompGetSubVecs(in,NULL,&inx);
382: VecCompGetSubVecs(out,NULL,&outx);
383: n = data->r_e - data->r_s;
384: SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
385: /* auxV <- in */
386: for (i=0;i<n;i++) {
387: VecCopy(inx[i],auxV[i]);
388: }
389: /* auxV <- auxV - u*(v'*auxV) */
390: dvd_improvex_applytrans_proj(data->d,auxV,n);
391: /* out <- K' * aux */
392: for (i=0;i<n;i++) {
393: PCApplyTranspose(data->old_pc,auxV[i],outx[i]);
394: }
395: SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
396: return(0);
397: }
399: static PetscErrorCode MatMult_dvd_jd(Mat A,Vec in,Vec out)
400: {
402: dvdImprovex_jd *data;
403: PetscInt n;
404: const Vec *inx, *outx;
405: PCSide side;
408: MatShellGetContext(A,(void**)&data);
409: VecCompGetSubVecs(in,NULL,&inx);
410: VecCompGetSubVecs(out,NULL,&outx);
411: n = data->r_e - data->r_s;
412: /* out <- theta[1]A*in - theta[0]*B*in */
413: dvd_aux_matmult(data,inx,outx);
414: KSPGetPCSide(data->ksp,&side);
415: if (side == PC_RIGHT) {
416: /* out <- out - v*(u'*out) */
417: dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
418: }
419: return(0);
420: }
422: static PetscErrorCode MatMultTranspose_dvd_jd(Mat A,Vec in,Vec out)
423: {
425: dvdImprovex_jd *data;
426: PetscInt n,i;
427: const Vec *inx,*outx,*r;
428: Vec *auxV;
429: PCSide side;
432: MatShellGetContext(A,(void**)&data);
433: VecCompGetSubVecs(in,NULL,&inx);
434: VecCompGetSubVecs(out,NULL,&outx);
435: n = data->r_e - data->r_s;
436: KSPGetPCSide(data->ksp,&side);
437: if (side == PC_RIGHT) {
438: /* auxV <- in */
439: SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
440: for (i=0;i<n;i++) {
441: VecCopy(inx[i],auxV[i]);
442: }
443: /* auxV <- auxV - v*(u'*auxV) */
444: dvd_improvex_applytrans_proj(data->d,auxV,n);
445: r = auxV;
446: } else r = inx;
447: /* out <- theta[1]A*r - theta[0]*B*r */
448: dvd_aux_matmulttrans(data,r,outx);
449: if (side == PC_RIGHT) {
450: SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
451: }
452: return(0);
453: }
455: static PetscErrorCode MatCreateVecs_dvd_jd(Mat A,Vec *right,Vec *left)
456: {
458: Vec *r,*l;
459: dvdImprovex_jd *data;
460: PetscInt n,i;
463: MatShellGetContext(A,(void**)&data);
464: n = data->ksp_max_size;
465: if (right) {
466: PetscMalloc1(n,&r);
467: }
468: if (left) {
469: PetscMalloc1(n,&l);
470: }
471: for (i=0;i<n;i++) {
472: MatCreateVecs(data->d->A,right?&r[i]:NULL,left?&l[i]:NULL);
473: }
474: if (right) {
475: VecCreateCompWithVecs(r,n,data->friends,right);
476: for (i=0;i<n;i++) {
477: VecDestroy(&r[i]);
478: }
479: }
480: if (left) {
481: VecCreateCompWithVecs(l,n,data->friends,left);
482: for (i=0;i<n;i++) {
483: VecDestroy(&l[i]);
484: }
485: }
487: if (right) {
488: PetscFree(r);
489: }
490: if (left) {
491: PetscFree(l);
492: }
493: return(0);
494: }
496: static PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)
497: {
499: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
500: PetscInt rA, cA, rlA, clA;
501: Mat A;
502: PetscBool t;
503: PC pc;
504: Vec v0[2];
507: data->size_cX = data->old_size_X = 0;
508: data->lastTol = data->dynamic?0.5:0.0;
510: /* Setup the ksp */
511: if (data->ksp) {
512: /* Create the reference vector */
513: BVGetColumn(d->eps->V,0,&v0[0]);
514: v0[1] = v0[0];
515: VecCreateCompWithVecs(v0,data->ksp_max_size,NULL,&data->friends);
516: BVRestoreColumn(d->eps->V,0,&v0[0]);
517: PetscLogObjectParent((PetscObject)d->eps,(PetscObject)data->friends);
519: /* Save the current pc and set a PCNONE */
520: KSPGetPC(data->ksp, &data->old_pc);
521: PetscObjectTypeCompare((PetscObject)data->old_pc,PCNONE,&t);
522: data->lastTol = 0.5;
523: if (t) data->old_pc = 0;
524: else {
525: PetscObjectReference((PetscObject)data->old_pc);
526: PCCreate(PetscObjectComm((PetscObject)d->eps),&pc);
527: PCSetType(pc,PCSHELL);
528: PCSetOperators(pc,d->A,d->A);
529: PCSetReusePreconditioner(pc,PETSC_TRUE);
530: PCShellSetApply(pc,PCApply_dvd);
531: PCShellSetApplyBA(pc,PCApplyBA_dvd);
532: PCShellSetApplyTranspose(pc,PCApplyTranspose_dvd);
533: KSPSetPC(data->ksp,pc);
534: PCDestroy(&pc);
535: }
537: /* Create the (I-v*u')*K*(A-s*B) matrix */
538: MatGetSize(d->A,&rA,&cA);
539: MatGetLocalSize(d->A,&rlA,&clA);
540: 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);
541: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_dvd_jd);
542: MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_dvd_jd);
543: MatShellSetOperation(A,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_dvd_jd);
545: /* Try to avoid KSPReset */
546: KSPGetOperatorsSet(data->ksp,&t,NULL);
547: if (t) {
548: Mat M;
549: PetscInt rM;
550: KSPGetOperators(data->ksp,&M,NULL);
551: MatGetSize(M,&rM,NULL);
552: if (rM != rA*data->ksp_max_size) { KSPReset(data->ksp); }
553: }
554: KSPSetOperators(data->ksp,A,A);
555: KSPSetReusePreconditioner(data->ksp,PETSC_TRUE);
556: KSPSetUp(data->ksp);
557: MatDestroy(&A);
558: } else {
559: data->old_pc = 0;
560: data->friends = NULL;
561: }
562: BVSetActiveColumns(data->KZ,0,0);
563: BVSetActiveColumns(data->U,0,0);
564: return(0);
565: }
567: /*
568: Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
569: kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
570: where
571: pX,pY, the right and left eigenvectors of the projected system
572: ld, the leading dimension of pX and pY
573: */
574: 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)
575: {
576: #if defined(PETSC_MISSING_LAPACK_GETRF)
578: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable");
579: #else
581: PetscInt n=i_e-i_s,size_KZ,V_new,rm,i,lv,kv,lKZ,kKZ;
582: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
583: PetscScalar *array;
584: Mat M;
585: Vec u[2],v[2];
586: PetscBLASInt s,ldXKZ,info;
589: /* Check consistency */
590: BVGetActiveColumns(d->eps->V,&lv,&kv);
591: V_new = lv - data->size_cX;
592: if (V_new > data->old_size_X) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
593: data->old_size_X = n;
594: data->size_cX = lv;
596: /* KZ <- KZ(rm:rm+max_cX-1) */
597: BVGetActiveColumns(data->KZ,&lKZ,&kKZ);
598: rm = PetscMax(V_new+lKZ,0);
599: if (rm > 0) {
600: for (i=0;i<lKZ;i++) {
601: BVCopyColumn(data->KZ,i+rm,i);
602: BVCopyColumn(data->U,i+rm,i);
603: }
604: }
606: /* XKZ <- XKZ(rm:rm+max_cX-1,rm:rm+max_cX-1) */
607: if (rm > 0) {
608: for (i=0;i<lKZ;i++) {
609: PetscMemcpy(&data->XKZ[i*data->ldXKZ+i],&data->XKZ[(i+rm)*data->ldXKZ+i+rm],sizeof(PetscScalar)*lKZ);
610: }
611: }
612: lKZ = PetscMin(0,lKZ+V_new);
613: BVSetActiveColumns(data->KZ,lKZ,lKZ+n);
614: BVSetActiveColumns(data->U,lKZ,lKZ+n);
616: /* Compute X, KZ and KR */
617: BVGetColumn(data->U,lKZ,u);
618: if (n>1) { BVGetColumn(data->U,lKZ+1,&u[1]); }
619: BVGetColumn(data->KZ,lKZ,v);
620: if (n>1) { BVGetColumn(data->KZ,lKZ+1,&v[1]); }
621: d->improvex_jd_proj_uv(d,i_s,i_e,u,v,kr,theta,thetai,pX,pY,ld);
622: BVRestoreColumn(data->U,lKZ,u);
623: if (n>1) { BVRestoreColumn(data->U,lKZ+1,&u[1]); }
624: BVRestoreColumn(data->KZ,lKZ,v);
625: if (n>1) { BVRestoreColumn(data->KZ,lKZ+1,&v[1]); }
627: /* XKZ <- U'*KZ */
628: MatCreateSeqDense(PETSC_COMM_SELF,lKZ+n,lKZ+n,NULL,&M);
629: BVMatProject(data->KZ,NULL,data->U,M);
630: MatDenseGetArray(M,&array);
631: for (i=lKZ;i<lKZ+n;i++) { /* upper part */
632: PetscMemcpy(&data->XKZ[data->ldXKZ*i],&array[i*(lKZ+n)],lKZ*sizeof(PetscScalar));
633: }
634: for (i=0;i<lKZ+n;i++) { /* lower part */
635: PetscMemcpy(&data->XKZ[data->ldXKZ*i+lKZ],&array[i*(lKZ+n)+lKZ],n*sizeof(PetscScalar));
636: }
637: MatDenseRestoreArray(M,&array);
638: MatDestroy(&M);
640: /* iXKZ <- inv(XKZ) */
641: size_KZ = lKZ+n;
642: PetscBLASIntCast(lKZ+n,&s);
643: data->ldiXKZ = data->size_iXKZ = size_KZ;
644: for (i=0;i<size_KZ;i++) {
645: PetscMemcpy(&data->iXKZ[data->ldiXKZ*i],&data->XKZ[data->ldXKZ*i],sizeof(PetscScalar)*size_KZ);
646: }
647: PetscBLASIntCast(data->ldiXKZ,&ldXKZ);
648: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
649: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&s,&s,data->iXKZ,&ldXKZ,data->iXKZPivots,&info));
650: PetscFPTrapPop();
651: SlepcCheckLapackInfo("getrf",info);
652: return(0);
653: #endif
654: }
656: static PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
657: {
658: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
660: PetscInt i,j,n,maxits,maxits0,lits,s,ld,k,max_size_D,lV,kV;
661: PetscScalar *pX,*pY;
662: PetscReal tol,tol0;
663: Vec *kr,kr_comp,D_comp,D[2],kr0[2];
664: PetscBool odd_situation = PETSC_FALSE;
667: BVGetActiveColumns(d->eps->V,&lV,&kV);
668: max_size_D = d->eps->ncv-kV;
669: /* Quick exit */
670: if ((max_size_D == 0) || r_e-r_s <= 0) {
671: *size_D = 0;
672: return(0);
673: }
675: n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);
676: if (n == 0) SETERRQ(PETSC_COMM_SELF,1,"n == 0");
677: if (data->size_X < r_e-r_s) SETERRQ(PETSC_COMM_SELF,1,"size_X < r_e-r_s");
679: DSGetLeadingDimension(d->eps->ds,&ld);
681: /* Restart lastTol if a new pair converged */
682: if (data->dynamic && data->size_cX < lV)
683: data->lastTol = 0.5;
685: for (i=0,s=0;i<n;i+=s) {
686: /* If the selected eigenvalue is complex, but the arithmetic is real... */
687: #if !defined(PETSC_USE_COMPLEX)
688: if (d->eigi[r_s+i] != 0.0) {
689: if (i+2 <= max_size_D) s=2;
690: else break;
691: } else
692: #endif
693: s=1;
695: data->r_s = r_s+i;
696: data->r_e = r_s+i+s;
697: SlepcVecPoolGetVecs(d->auxV,s,&kr);
699: /* Compute theta, maximum iterations and tolerance */
700: maxits = 0;
701: tol = 1;
702: for (j=0;j<s;j++) {
703: d->improvex_jd_lit(d,r_s+i+j,&data->theta[2*j],&data->thetai[j],&maxits0,&tol0);
704: maxits += maxits0;
705: tol *= tol0;
706: }
707: maxits/= s;
708: tol = data->dynamic?data->lastTol:PetscExpReal(PetscLogReal(tol)/s);
710: /* Compute u, v and kr */
711: k = r_s+i;
712: DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);
713: k = r_s+i;
714: DSVectors(d->eps->ds,DS_MAT_Y,&k,NULL);
715: DSGetArray(d->eps->ds,DS_MAT_X,&pX);
716: DSGetArray(d->eps->ds,DS_MAT_Y,&pY);
717: dvd_improvex_jd_proj_cuv(d,r_s+i,r_s+i+s,kr,data->theta,data->thetai,pX,pY,ld);
718: DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);
719: DSRestoreArray(d->eps->ds,DS_MAT_Y,&pY);
721: /* Check if the first eigenpairs are converged */
722: if (i == 0) {
723: PetscInt oldnpreconv = d->npreconv;
724: d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv);
725: if (d->npreconv > oldnpreconv) break;
726: }
728: /* Test the odd situation of solving Ax=b with A=I */
729: #if !defined(PETSC_USE_COMPLEX)
730: odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && data->thetai[0] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
731: #else
732: odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
733: #endif
734: /* If JD */
735: if (data->ksp && !odd_situation) {
736: /* kr <- -kr */
737: for (j=0;j<s;j++) {
738: VecScale(kr[j],-1.0);
739: }
741: /* Compose kr and D */
742: kr0[0] = kr[0];
743: kr0[1] = (s==2 ? kr[1] : NULL);
744: VecCreateCompWithVecs(kr0,data->ksp_max_size,data->friends,&kr_comp);
745: BVGetColumn(d->eps->V,kV+i,&D[0]);
746: if (s==2) { BVGetColumn(d->eps->V,kV+i+1,&D[1]); }
747: else D[1] = NULL;
748: VecCreateCompWithVecs(D,data->ksp_max_size,data->friends,&D_comp);
749: VecCompSetSubVecs(data->friends,s,NULL);
751: /* Solve the correction equation */
752: KSPSetTolerances(data->ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,maxits);
753: KSPSolve(data->ksp,kr_comp,D_comp);
754: KSPGetIterationNumber(data->ksp,&lits);
756: /* Destroy the composed ks and D */
757: VecDestroy(&kr_comp);
758: VecDestroy(&D_comp);
759: BVRestoreColumn(d->eps->V,kV+i,&D[0]);
760: if (s==2) { BVRestoreColumn(d->eps->V,kV+i+1,&D[1]); }
762: /* If GD */
763: } else {
764: BVGetColumn(d->eps->V,kV+i,&D[0]);
765: if (s==2) { BVGetColumn(d->eps->V,kV+i+1,&D[1]); }
766: for (j=0;j<s;j++) {
767: d->improvex_precond(d,r_s+i+j,kr[j],D[j]);
768: }
769: dvd_improvex_apply_proj(d,D,s);
770: BVRestoreColumn(d->eps->V,kV+i,&D[0]);
771: if (s==2) { BVRestoreColumn(d->eps->V,kV+i+1,&D[1]); }
772: }
773: /* Prevent that short vectors are discarded in the orthogonalization */
774: if (i == 0 && d->eps->errest[d->nconv+r_s] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s] < PETSC_MAX_REAL) {
775: for (j=0;j<s;j++) {
776: BVScaleColumn(d->eps->V,kV+i+j,1.0/d->eps->errest[d->nconv+r_s]);
777: }
778: }
779: SlepcVecPoolRestoreVecs(d->auxV,s,&kr);
780: }
781: *size_D = i;
782: if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);
783: return(0);
784: }
786: PetscErrorCode dvd_improvex_jd(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs,PetscBool dynamic)
787: {
789: dvdImprovex_jd *data;
790: PetscBool useGD;
791: PC pc;
792: PetscInt size_P;
795: /* Setting configuration constrains */
796: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&useGD);
798: /* If the arithmetic is real and the problem is not Hermitian, then
799: the block size is incremented in one */
800: #if !defined(PETSC_USE_COMPLEX)
801: if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) {
802: max_bs++;
803: b->max_size_P = PetscMax(b->max_size_P,2);
804: } else
805: #endif
806: {
807: b->max_size_P = PetscMax(b->max_size_P,1);
808: }
809: b->max_size_X = PetscMax(b->max_size_X,max_bs);
810: size_P = b->max_size_P;
812: /* Setup the preconditioner */
813: if (ksp) {
814: KSPGetPC(ksp,&pc);
815: dvd_static_precond_PC(d,b,pc);
816: } else {
817: dvd_static_precond_PC(d,b,0);
818: }
820: /* Setup the step */
821: if (b->state >= DVD_STATE_CONF) {
822: PetscNewLog(d->eps,&data);
823: data->dynamic = dynamic;
824: PetscMalloc1(size_P*size_P,&data->XKZ);
825: PetscMalloc1(size_P*size_P,&data->iXKZ);
826: PetscMalloc1(size_P,&data->iXKZPivots);
827: data->ldXKZ = size_P;
828: data->size_X = b->max_size_X;
829: d->improveX_data = data;
830: data->ksp = useGD? NULL: ksp;
831: data->d = d;
832: d->improveX = dvd_improvex_jd_gen;
833: #if !defined(PETSC_USE_COMPLEX)
834: if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) data->ksp_max_size = 2;
835: else
836: #endif
837: data->ksp_max_size = 1;
838: /* Create various vector basis */
839: BVDuplicateResize(d->eps->V,size_P,&data->KZ);
840: BVSetMatrix(data->KZ,NULL,PETSC_FALSE);
841: BVDuplicate(data->KZ,&data->U);
843: EPSDavidsonFLAdd(&d->startList,dvd_improvex_jd_start);
844: EPSDavidsonFLAdd(&d->endList,dvd_improvex_jd_end);
845: EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_jd_d);
846: }
847: return(0);
848: }
850: #if !defined(PETSC_USE_COMPLEX)
851: PETSC_STATIC_INLINE PetscErrorCode dvd_complex_rayleigh_quotient(Vec ur,Vec ui,Vec Axr,Vec Axi,Vec Bxr,Vec Bxi,PetscScalar *eigr,PetscScalar *eigi)
852: {
854: PetscScalar rAr,iAr,rAi,iAi,rBr,iBr,rBi,iBi,b0,b2,b4,b6,b7;
857: /* eigr = [(rAr+iAi)*(rBr+iBi) + (rAi-iAr)*(rBi-iBr)]/k
858: eigi = [(rAi-iAr)*(rBr+iBi) - (rAr+iAi)*(rBi-iBr)]/k
859: k = (rBr+iBi)*(rBr+iBi) + (rBi-iBr)*(rBi-iBr) */
860: VecDotBegin(Axr,ur,&rAr); /* r*A*r */
861: VecDotBegin(Axr,ui,&iAr); /* i*A*r */
862: VecDotBegin(Axi,ur,&rAi); /* r*A*i */
863: VecDotBegin(Axi,ui,&iAi); /* i*A*i */
864: VecDotBegin(Bxr,ur,&rBr); /* r*B*r */
865: VecDotBegin(Bxr,ui,&iBr); /* i*B*r */
866: VecDotBegin(Bxi,ur,&rBi); /* r*B*i */
867: VecDotBegin(Bxi,ui,&iBi); /* i*B*i */
868: VecDotEnd(Axr,ur,&rAr); /* r*A*r */
869: VecDotEnd(Axr,ui,&iAr); /* i*A*r */
870: VecDotEnd(Axi,ur,&rAi); /* r*A*i */
871: VecDotEnd(Axi,ui,&iAi); /* i*A*i */
872: VecDotEnd(Bxr,ur,&rBr); /* r*B*r */
873: VecDotEnd(Bxr,ui,&iBr); /* i*B*r */
874: VecDotEnd(Bxi,ur,&rBi); /* r*B*i */
875: VecDotEnd(Bxi,ui,&iBi); /* i*B*i */
876: b0 = rAr+iAi; /* rAr+iAi */
877: b2 = rAi-iAr; /* rAi-iAr */
878: b4 = rBr+iBi; /* rBr+iBi */
879: b6 = rBi-iBr; /* rBi-iBr */
880: b7 = b4*b4 + b6*b6; /* k */
881: *eigr = (b0*b4 + b2*b6) / b7; /* eig_r */
882: *eigi = (b2*b4 - b0*b6) / b7; /* eig_i */
883: return(0);
884: }
885: #endif
887: PETSC_STATIC_INLINE PetscErrorCode dvd_compute_n_rr(PetscInt i_s,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,Vec *u,Vec *Ax,Vec *Bx)
888: {
890: PetscInt i;
891: PetscScalar b0,b1;
894: for (i=0; i<n; i++) {
895: #if !defined(PETSC_USE_COMPLEX)
896: if (eigi[i_s+i] != 0.0) {
897: PetscScalar eigr0=0.0,eigi0=0.0;
898: dvd_complex_rayleigh_quotient(u[i],u[i+1],Ax[i],Ax[i+1],Bx[i],Bx[i+1],&eigr0,&eigi0);
899: 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) {
900: PetscInfo4(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);
901: }
902: i++;
903: } else
904: #endif
905: {
906: VecDotBegin(Ax[i],u[i],&b0);
907: VecDotBegin(Bx[i],u[i],&b1);
908: VecDotEnd(Ax[i],u[i],&b0);
909: VecDotEnd(Bx[i],u[i],&b1);
910: b0 = b0/b1;
911: if (PetscAbsScalar(eigr[i_s+i]-b0)/PetscAbsScalar(eigr[i_s+i]) > 1e-10) {
912: PetscInfo4(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));
913: }
914: }
915: }
916: return(0);
917: }
919: /*
920: Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
921: kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
922: where
923: pX,pY, the right and left eigenvectors of the projected system
924: ld, the leading dimension of pX and pY
925: */
926: 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)
927: {
929: PetscInt n = i_e-i_s,i;
930: PetscScalar *b;
931: Vec *Ax,*Bx,*r;
932: Mat M;
933: BV X;
936: BVDuplicateResize(d->eps->V,4,&X);
937: MatCreateSeqDense(PETSC_COMM_SELF,4,4,NULL,&M);
938: /* u <- X(i) */
939: dvd_improvex_compute_X(d,i_s,i_e,u,pX,ld);
941: /* v <- theta[0]A*u + theta[1]*B*u */
943: /* Bx <- B*X(i) */
944: Bx = kr;
945: if (d->BX) {
946: for (i=i_s; i<i_e; ++i) {
947: BVMultVec(d->BX,1.0,0.0,Bx[i-i_s],&pX[ld*i]);
948: }
949: } else {
950: for (i=0;i<n;i++) {
951: if (d->B) {
952: MatMult(d->B, u[i], Bx[i]);
953: } else {
954: VecCopy(u[i], Bx[i]);
955: }
956: }
957: }
959: /* Ax <- A*X(i) */
960: SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&r);
961: Ax = r;
962: for (i=i_s; i<i_e; ++i) {
963: BVMultVec(d->AX,1.0,0.0,Ax[i-i_s],&pX[ld*i]);
964: }
966: /* v <- Y(i) */
967: for (i=i_s; i<i_e; ++i) {
968: BVMultVec(d->W?d->W:d->eps->V,1.0,0.0,v[i-i_s],&pY[ld*i]);
969: }
971: /* Recompute the eigenvalue */
972: dvd_compute_n_rr(i_s,n,d->eigr,d->eigi,v,Ax,Bx);
974: for (i=0;i<n;i++) {
975: #if !defined(PETSC_USE_COMPLEX)
976: if (d->eigi[i_s+i] != 0.0) {
977: /* [r_i r_i+1 kr_i kr_i+1]*= [ theta_2i' 0 1 0
978: 0 theta_2i' 0 1
979: theta_2i+1 -thetai_i -eigr_i -eigi_i
980: thetai_i theta_2i+1 eigi_i -eigr_i ] */
981: MatDenseGetArray(M,&b);
982: b[0] = b[5] = PetscConj(theta[2*i]);
983: b[2] = b[7] = -theta[2*i+1];
984: b[6] = -(b[3] = thetai[i]);
985: b[1] = b[4] = 0.0;
986: b[8] = b[13] = 1.0/d->nX[i_s+i];
987: b[10] = b[15] = -d->eigr[i_s+i]/d->nX[i_s+i];
988: b[14] = -(b[11] = d->eigi[i_s+i]/d->nX[i_s+i]);
989: b[9] = b[12] = 0.0;
990: MatDenseRestoreArray(M,&b);
991: BVInsertVec(X,0,Ax[i]);
992: BVInsertVec(X,1,Ax[i+1]);
993: BVInsertVec(X,2,Bx[i]);
994: BVInsertVec(X,3,Bx[i+1]);
995: BVSetActiveColumns(X,0,4);
996: BVMultInPlace(X,M,0,4);
997: BVCopyVec(X,0,Ax[i]);
998: BVCopyVec(X,1,Ax[i+1]);
999: BVCopyVec(X,2,Bx[i]);
1000: BVCopyVec(X,3,Bx[i+1]);
1001: i++;
1002: } else
1003: #endif
1004: {
1005: /* [Ax_i Bx_i]*= [ theta_2i' 1/nX_i
1006: theta_2i+1 -eig_i/nX_i ] */
1007: MatDenseGetArray(M,&b);
1008: b[0] = PetscConj(theta[i*2]);
1009: b[1] = theta[i*2+1];
1010: b[4] = 1.0/d->nX[i_s+i];
1011: b[5] = -d->eigr[i_s+i]/d->nX[i_s+i];
1012: MatDenseRestoreArray(M,&b);
1013: BVInsertVec(X,0,Ax[i]);
1014: BVInsertVec(X,1,Bx[i]);
1015: BVSetActiveColumns(X,0,2);
1016: BVMultInPlace(X,M,0,2);
1017: BVCopyVec(X,0,Ax[i]);
1018: BVCopyVec(X,1,Bx[i]);
1019: }
1020: }
1021: for (i=0; i<n; i++) d->nX[i_s+i] = 1.0;
1023: /* v <- K^{-1} r = K^{-1}(theta_2i'*Ax + theta_2i+1*Bx) */
1024: for (i=0;i<n;i++) {
1025: d->improvex_precond(d,i_s+i,r[i],v[i]);
1026: }
1027: SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&r);
1029: /* kr <- P*(Ax - eig_i*Bx) */
1030: d->calcpairs_proj_res(d,i_s,i_e,kr);
1031: BVDestroy(&X);
1032: MatDestroy(&M);
1033: return(0);
1034: }
1036: static PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d,PetscInt i,PetscScalar* theta,PetscScalar* thetai,PetscInt *maxits,PetscReal *tol)
1037: {
1038: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
1039: PetscReal a;
1042: a = SlepcAbsEigenvalue(d->eigr[i],d->eigi[i]);
1044: if (d->nR[i] < data->fix*a) {
1045: theta[0] = d->eigr[i];
1046: theta[1] = 1.0;
1047: #if !defined(PETSC_USE_COMPLEX)
1048: *thetai = d->eigi[i];
1049: #endif
1050: } else {
1051: theta[0] = d->target[0];
1052: theta[1] = d->target[1];
1053: #if !defined(PETSC_USE_COMPLEX)
1054: *thetai = 0.0;
1055: #endif
1056: }
1058: #if defined(PETSC_USE_COMPLEX)
1059: if (thetai) *thetai = 0.0;
1060: #endif
1061: *maxits = data->maxits;
1062: *tol = data->tol;
1063: return(0);
1064: }
1066: PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard *d,dvdBlackboard *b,PetscInt maxits,PetscReal tol,PetscReal fix)
1067: {
1068: dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
1071: /* Setup the step */
1072: if (b->state >= DVD_STATE_CONF) {
1073: data->maxits = maxits;
1074: data->tol = tol;
1075: data->fix = fix;
1076: d->improvex_jd_lit = dvd_improvex_jd_lit_const_0;
1077: }
1078: return(0);
1079: }
1081: PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d,dvdBlackboard *b)
1082: {
1084: /* Setup the step */
1085: if (b->state >= DVD_STATE_CONF) {
1086: d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX;
1087: }
1088: return(0);
1089: }
1091: PetscErrorCode dvd_improvex_compute_X(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u_,PetscScalar *pX,PetscInt ld)
1092: {
1094: PetscInt n = i_e - i_s,i;
1095: Vec *u;
1098: if (u_) u = u_;
1099: else if (d->correctXnorm) {
1100: SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&u);
1101: }
1102: if (u_ || d->correctXnorm) {
1103: for (i=0; i<n; i++) {
1104: BVMultVec(d->eps->V,1.0,0.0,u[i],&pX[ld*(i+i_s)]);
1105: }
1106: }
1107: /* nX(i) <- ||X(i)|| */
1108: if (d->correctXnorm) {
1109: for (i=0; i<n; i++) {
1110: VecNormBegin(u[i],NORM_2,&d->nX[i_s+i]);
1111: }
1112: for (i=0; i<n; i++) {
1113: VecNormEnd(u[i],NORM_2,&d->nX[i_s+i]);
1114: }
1115: #if !defined(PETSC_USE_COMPLEX)
1116: for (i=0;i<n;i++) {
1117: if (d->eigi[i_s+i] != 0.0) {
1118: 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]);
1119: i++;
1120: }
1121: }
1122: #endif
1123: } else {
1124: for (i=0;i<n;i++) d->nX[i_s+i] = 1.0;
1125: }
1126: if (d->correctXnorm && !u_) {
1127: SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&u);
1128: }
1129: return(0);
1130: }