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