Actual source code: dvdimprovex.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  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:   PetscAssert(cV<=2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");

 64:   /* h <- X'*V */
 65:   PetscMalloc1(data->size_iXKZ*cV,&h);
 66:   ldh = data->size_iXKZ;
 67:   BVGetActiveColumns(data->U,&l,&k);
 68:   PetscAssert(ldh==k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
 69:   BVSetActiveColumns(data->U,0,k);
 70:   for (i=0;i<cV;i++) {
 71:     BVDotVec(data->U,V[i],&h[ldh*i]);
 72: #if defined(PETSC_USE_COMPLEX)
 73:     for (j=0; j<k; j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
 74: #endif
 75:   }
 76:   BVSetActiveColumns(data->U,l,k);

 78:   /* h <- iXKZ\h */
 79:   PetscBLASIntCast(cV,&cV_);
 80:   PetscBLASIntCast(data->size_iXKZ,&n);
 81:   PetscBLASIntCast(data->ldiXKZ,&ld);
 82:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 83:   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
 84:   PetscFPTrapPop();
 85:   SlepcCheckLapackInfo("getrs",info);

 87:   /* V <- V - KZ*h */
 88:   BVSetActiveColumns(data->KZ,0,k);
 89:   for (i=0;i<cV;i++) BVMultVec(data->KZ,-1.0,1.0,V[i],&h[ldh*i]);
 90:   BVSetActiveColumns(data->KZ,l,k);
 91:   PetscFree(h);
 92:   PetscFunctionReturn(0);
 93: }

 95: /*
 96:    Compute (I - X*iXKZ*KZ')*V where,
 97:    V, the vectors to apply the projector,
 98:    cV, the number of vectors in V,
 99: */
100: static PetscErrorCode dvd_improvex_applytrans_proj(dvdDashboard *d,Vec *V,PetscInt cV)
101: {
102:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
103:   PetscInt       i,ldh,k,l;
104:   PetscScalar    *h;
105:   PetscBLASInt   cV_, n, info, ld;
106: #if defined(PETSC_USE_COMPLEX)
107:   PetscInt       j;
108: #endif

110:   PetscAssert(cV<=2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");

112:   /* h <- KZ'*V */
113:   PetscMalloc1(data->size_iXKZ*cV,&h);
114:   ldh = data->size_iXKZ;
115:   BVGetActiveColumns(data->U,&l,&k);
116:   PetscAssert(ldh==k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
117:   BVSetActiveColumns(data->KZ,0,k);
118:   for (i=0;i<cV;i++) {
119:     BVDotVec(data->KZ,V[i],&h[ldh*i]);
120: #if defined(PETSC_USE_COMPLEX)
121:     for (j=0;j<k;j++) h[ldh*i+j] = PetscConj(h[ldh*i+j]);
122: #endif
123:   }
124:   BVSetActiveColumns(data->KZ,l,k);

126:   /* h <- iXKZ\h */
127:   PetscBLASIntCast(cV,&cV_);
128:   PetscBLASIntCast(data->size_iXKZ,&n);
129:   PetscBLASIntCast(data->ldiXKZ,&ld);
130:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
131:   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
132:   PetscFPTrapPop();
133:   SlepcCheckLapackInfo("getrs",info);

135:   /* V <- V - U*h */
136:   BVSetActiveColumns(data->U,0,k);
137:   for (i=0;i<cV;i++) BVMultVec(data->U,-1.0,1.0,V[i],&h[ldh*i]);
138:   BVSetActiveColumns(data->U,l,k);
139:   PetscFree(h);
140:   PetscFunctionReturn(0);
141: }

143: static PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)
144: {
145:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;

147:   VecDestroy(&data->friends);

149:   /* Restore the pc of ksp */
150:   if (data->old_pc) {
151:     KSPSetPC(data->ksp, data->old_pc);
152:     PCDestroy(&data->old_pc);
153:   }
154:   PetscFunctionReturn(0);
155: }

157: static PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)
158: {
159:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;

161:   /* Free local data and objects */
162:   PetscFree(data->XKZ);
163:   PetscFree(data->iXKZ);
164:   PetscFree(data->iXKZPivots);
165:   BVDestroy(&data->KZ);
166:   BVDestroy(&data->U);
167:   PetscFree(data);
168:   PetscFunctionReturn(0);
169: }

171: /*
172:    y <- theta[1]A*x - theta[0]*B*x
173:    auxV, two auxiliary vectors
174:  */
175: static inline PetscErrorCode dvd_aux_matmult(dvdImprovex_jd *data,const Vec *x,const Vec *y)
176: {
177:   PetscInt       n,i;
178:   const Vec      *Bx;
179:   Vec            *auxV;

181:   n = data->r_e - data->r_s;
182:   for (i=0;i<n;i++) MatMult(data->d->A,x[i],y[i]);

184:   SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);
185:   for (i=0;i<n;i++) {
186: #if !defined(PETSC_USE_COMPLEX)
187:     if (PetscUnlikely(data->d->eigi[data->r_s+i] != 0.0)) {
188:       if (data->d->B) {
189:         MatMult(data->d->B,x[i],auxV[0]);
190:         MatMult(data->d->B,x[i+1],auxV[1]);
191:         Bx = auxV;
192:       } else Bx = &x[i];

194:       /* y_i   <- [ t_2i+1*A*x_i   - t_2i*Bx_i + ti_i*Bx_i+1;
195:          y_i+1      t_2i+1*A*x_i+1 - ti_i*Bx_i - t_2i*Bx_i+1  ] */
196:       VecAXPBYPCZ(y[i],-data->theta[2*i],data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);
197:       VecAXPBYPCZ(y[i+1],-data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);
198:       i++;
199:     } else
200: #endif
201:     {
202:       if (data->d->B) {
203:         MatMult(data->d->B,x[i],auxV[0]);
204:         Bx = auxV;
205:       } else Bx = &x[i];
206:       VecAXPBY(y[i],-data->theta[i*2],data->theta[i*2+1],Bx[0]);
207:     }
208:   }
209:   SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);
210:   PetscFunctionReturn(0);
211: }

213: /*
214:    y <- theta[1]'*A'*x - theta[0]'*B'*x
215:  */
216: static inline PetscErrorCode dvd_aux_matmulttrans(dvdImprovex_jd *data,const Vec *x,const Vec *y)
217: {
218:   PetscInt       n,i;
219:   const Vec      *Bx;
220:   Vec            *auxV;

222:   n = data->r_e - data->r_s;
223:   for (i=0;i<n;i++) MatMultTranspose(data->d->A,x[i],y[i]);

225:   SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);
226:   for (i=0;i<n;i++) {
227: #if !defined(PETSC_USE_COMPLEX)
228:     if (data->d->eigi[data->r_s+i] != 0.0) {
229:       if (data->d->B) {
230:         MatMultTranspose(data->d->B,x[i],auxV[0]);
231:         MatMultTranspose(data->d->B,x[i+1],auxV[1]);
232:         Bx = auxV;
233:       } else Bx = &x[i];

235:       /* y_i   <- [ t_2i+1*A*x_i   - t_2i*Bx_i - ti_i*Bx_i+1;
236:          y_i+1      t_2i+1*A*x_i+1 + ti_i*Bx_i - t_2i*Bx_i+1  ] */
237:       VecAXPBYPCZ(y[i],-data->theta[2*i],-data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);
238:       VecAXPBYPCZ(y[i+1],data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);
239:       i++;
240:     } else
241: #endif
242:     {
243:       if (data->d->B) {
244:         MatMultTranspose(data->d->B,x[i],auxV[0]);
245:         Bx = auxV;
246:       } else Bx = &x[i];
247:       VecAXPBY(y[i],PetscConj(-data->theta[i*2]),PetscConj(data->theta[i*2+1]),Bx[0]);
248:     }
249:   }
250:   SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);
251:   PetscFunctionReturn(0);
252: }

254: static PetscErrorCode PCApplyBA_dvd(PC pc,PCSide side,Vec in,Vec out,Vec w)
255: {
256:   dvdImprovex_jd *data;
257:   PetscInt       n,i;
258:   const Vec      *inx,*outx,*wx;
259:   Vec            *auxV;
260:   Mat            A;

262:   PCGetOperators(pc,&A,NULL);
263:   MatShellGetContext(A,&data);
264:   VecCompGetSubVecs(in,NULL,&inx);
265:   VecCompGetSubVecs(out,NULL,&outx);
266:   VecCompGetSubVecs(w,NULL,&wx);
267:   n = data->r_e - data->r_s;
268:   SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
269:   switch (side) {
270:   case PC_LEFT:
271:     /* aux <- theta[1]A*in - theta[0]*B*in */
272:     dvd_aux_matmult(data,inx,auxV);

274:     /* out <- K * aux */
275:     for (i=0;i<n;i++) data->d->improvex_precond(data->d,data->r_s+i,auxV[i],outx[i]);
276:     break;
277:   case PC_RIGHT:
278:     /* aux <- K * in */
279:     for (i=0;i<n;i++) data->d->improvex_precond(data->d,data->r_s+i,inx[i],auxV[i]);

281:     /* out <- theta[1]A*auxV - theta[0]*B*auxV */
282:     dvd_aux_matmult(data,auxV,outx);
283:     break;
284:   case PC_SYMMETRIC:
285:     /* aux <- K^{1/2} * in */
286:     for (i=0;i<n;i++) PCApplySymmetricRight(data->old_pc,inx[i],auxV[i]);

288:     /* wx <- theta[1]A*auxV - theta[0]*B*auxV */
289:     dvd_aux_matmult(data,auxV,wx);

291:     /* aux <- K^{1/2} * in */
292:     for (i=0;i<n;i++) PCApplySymmetricLeft(data->old_pc,wx[i],outx[i]);
293:     break;
294:   default:
295:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported KSP side");
296:   }
297:   /* out <- out - v*(u'*out) */
298:   dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
299:   SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
300:   PetscFunctionReturn(0);
301: }

303: static PetscErrorCode PCApply_dvd(PC pc,Vec in,Vec out)
304: {
305:   dvdImprovex_jd *data;
306:   PetscInt       n,i;
307:   const Vec      *inx, *outx;
308:   Mat            A;

310:   PCGetOperators(pc,&A,NULL);
311:   MatShellGetContext(A,&data);
312:   VecCompGetSubVecs(in,NULL,&inx);
313:   VecCompGetSubVecs(out,NULL,&outx);
314:   n = data->r_e - data->r_s;
315:   /* out <- K * in */
316:   for (i=0;i<n;i++) data->d->improvex_precond(data->d,data->r_s+i,inx[i],outx[i]);
317:   /* out <- out - v*(u'*out) */
318:   dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
319:   PetscFunctionReturn(0);
320: }

322: static PetscErrorCode PCApplyTranspose_dvd(PC pc,Vec in,Vec out)
323: {
324:   dvdImprovex_jd *data;
325:   PetscInt       n,i;
326:   const Vec      *inx, *outx;
327:   Vec            *auxV;
328:   Mat            A;

330:   PCGetOperators(pc,&A,NULL);
331:   MatShellGetContext(A,&data);
332:   VecCompGetSubVecs(in,NULL,&inx);
333:   VecCompGetSubVecs(out,NULL,&outx);
334:   n = data->r_e - data->r_s;
335:   SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
336:   /* auxV <- in */
337:   for (i=0;i<n;i++) VecCopy(inx[i],auxV[i]);
338:   /* auxV <- auxV - u*(v'*auxV) */
339:   dvd_improvex_applytrans_proj(data->d,auxV,n);
340:   /* out <- K' * aux */
341:   for (i=0;i<n;i++) PCApplyTranspose(data->old_pc,auxV[i],outx[i]);
342:   SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
343:   PetscFunctionReturn(0);
344: }

346: static PetscErrorCode MatMult_dvd_jd(Mat A,Vec in,Vec out)
347: {
348:   dvdImprovex_jd *data;
349:   PetscInt       n;
350:   const Vec      *inx, *outx;
351:   PCSide         side;

353:   MatShellGetContext(A,&data);
354:   VecCompGetSubVecs(in,NULL,&inx);
355:   VecCompGetSubVecs(out,NULL,&outx);
356:   n = data->r_e - data->r_s;
357:   /* out <- theta[1]A*in - theta[0]*B*in */
358:   dvd_aux_matmult(data,inx,outx);
359:   KSPGetPCSide(data->ksp,&side);
360:   if (side == PC_RIGHT) {
361:     /* out <- out - v*(u'*out) */
362:     dvd_improvex_apply_proj(data->d,(Vec*)outx,n);
363:   }
364:   PetscFunctionReturn(0);
365: }

367: static PetscErrorCode MatMultTranspose_dvd_jd(Mat A,Vec in,Vec out)
368: {
369:   dvdImprovex_jd *data;
370:   PetscInt       n,i;
371:   const Vec      *inx,*outx,*r;
372:   Vec            *auxV;
373:   PCSide         side;

375:   MatShellGetContext(A,&data);
376:   VecCompGetSubVecs(in,NULL,&inx);
377:   VecCompGetSubVecs(out,NULL,&outx);
378:   n = data->r_e - data->r_s;
379:   KSPGetPCSide(data->ksp,&side);
380:   if (side == PC_RIGHT) {
381:     /* auxV <- in */
382:     SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);
383:     for (i=0;i<n;i++) VecCopy(inx[i],auxV[i]);
384:     /* auxV <- auxV - v*(u'*auxV) */
385:     dvd_improvex_applytrans_proj(data->d,auxV,n);
386:     r = auxV;
387:   } else r = inx;
388:   /* out <- theta[1]A*r - theta[0]*B*r */
389:   dvd_aux_matmulttrans(data,r,outx);
390:   if (side == PC_RIGHT) SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);
391:   PetscFunctionReturn(0);
392: }

394: static PetscErrorCode MatCreateVecs_dvd_jd(Mat A,Vec *right,Vec *left)
395: {
396:   Vec            *r,*l;
397:   dvdImprovex_jd *data;
398:   PetscInt       n,i;

400:   MatShellGetContext(A,&data);
401:   n = data->ksp_max_size;
402:   if (right) PetscMalloc1(n,&r);
403:   if (left) PetscMalloc1(n,&l);
404:   for (i=0;i<n;i++) MatCreateVecs(data->d->A,right?&r[i]:NULL,left?&l[i]:NULL);
405:   if (right) {
406:     VecCreateCompWithVecs(r,n,data->friends,right);
407:     for (i=0;i<n;i++) VecDestroy(&r[i]);
408:   }
409:   if (left) {
410:     VecCreateCompWithVecs(l,n,data->friends,left);
411:     for (i=0;i<n;i++) VecDestroy(&l[i]);
412:   }

414:   if (right) PetscFree(r);
415:   if (left) PetscFree(l);
416:   PetscFunctionReturn(0);
417: }

419: static PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)
420: {
421:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
422:   PetscInt       rA, cA, rlA, clA;
423:   Mat            A;
424:   PetscBool      t;
425:   PC             pc;
426:   Vec            v0[2];

428:   data->size_cX = data->old_size_X = 0;
429:   data->lastTol = data->dynamic?0.5:0.0;

431:   /* Setup the ksp */
432:   if (data->ksp) {
433:     /* Create the reference vector */
434:     BVGetColumn(d->eps->V,0,&v0[0]);
435:     v0[1] = v0[0];
436:     VecCreateCompWithVecs(v0,data->ksp_max_size,NULL,&data->friends);
437:     BVRestoreColumn(d->eps->V,0,&v0[0]);
438:     PetscLogObjectParent((PetscObject)d->eps,(PetscObject)data->friends);

440:     /* Save the current pc and set a PCNONE */
441:     KSPGetPC(data->ksp, &data->old_pc);
442:     PetscObjectTypeCompare((PetscObject)data->old_pc,PCNONE,&t);
443:     data->lastTol = 0.5;
444:     if (t) data->old_pc = 0;
445:     else {
446:       PetscObjectReference((PetscObject)data->old_pc);
447:       PCCreate(PetscObjectComm((PetscObject)d->eps),&pc);
448:       PCSetType(pc,PCSHELL);
449:       PCSetOperators(pc,d->A,d->A);
450:       PCSetReusePreconditioner(pc,PETSC_TRUE);
451:       PCShellSetApply(pc,PCApply_dvd);
452:       PCShellSetApplyBA(pc,PCApplyBA_dvd);
453:       PCShellSetApplyTranspose(pc,PCApplyTranspose_dvd);
454:       KSPSetPC(data->ksp,pc);
455:       PCDestroy(&pc);
456:     }

458:     /* Create the (I-v*u')*K*(A-s*B) matrix */
459:     MatGetSize(d->A,&rA,&cA);
460:     MatGetLocalSize(d->A,&rlA,&clA);
461:     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);
462:     MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_dvd_jd);
463:     MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_dvd_jd);
464:     MatShellSetOperation(A,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_dvd_jd);

466:     /* Try to avoid KSPReset */
467:     KSPGetOperatorsSet(data->ksp,&t,NULL);
468:     if (t) {
469:       Mat      M;
470:       PetscInt rM;
471:       KSPGetOperators(data->ksp,&M,NULL);
472:       MatGetSize(M,&rM,NULL);
473:       if (rM != rA*data->ksp_max_size) KSPReset(data->ksp);
474:     }
475:     EPS_KSPSetOperators(data->ksp,A,A);
476:     KSPSetReusePreconditioner(data->ksp,PETSC_TRUE);
477:     KSPSetUp(data->ksp);
478:     MatDestroy(&A);
479:   } else {
480:     data->old_pc = 0;
481:     data->friends = NULL;
482:   }
483:   BVSetActiveColumns(data->KZ,0,0);
484:   BVSetActiveColumns(data->U,0,0);
485:   PetscFunctionReturn(0);
486: }

488: /*
489:   Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
490:   kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
491:   where
492:   pX,pY, the right and left eigenvectors of the projected system
493:   ld, the leading dimension of pX and pY
494: */
495: 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)
496: {
497:   PetscInt          n=i_e-i_s,size_KZ,V_new,rm,i,lv,kv,lKZ,kKZ;
498:   dvdImprovex_jd    *data = (dvdImprovex_jd*)d->improveX_data;
499:   const PetscScalar *array;
500:   Mat               M;
501:   Vec               u[2],v[2];
502:   PetscBLASInt      s,ldXKZ,info;

504:   /* Check consistency */
505:   BVGetActiveColumns(d->eps->V,&lv,&kv);
506:   V_new = lv - data->size_cX;
507:   PetscAssert(V_new<=data->old_size_X,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
508:   data->old_size_X = n;
509:   data->size_cX = lv;

511:   /* KZ <- KZ(rm:rm+max_cX-1) */
512:   BVGetActiveColumns(data->KZ,&lKZ,&kKZ);
513:   rm = PetscMax(V_new+lKZ,0);
514:   if (rm > 0) {
515:     for (i=0;i<lKZ;i++) {
516:       BVCopyColumn(data->KZ,i+rm,i);
517:       BVCopyColumn(data->U,i+rm,i);
518:     }
519:   }

521:   /* XKZ <- XKZ(rm:rm+max_cX-1,rm:rm+max_cX-1) */
522:   if (rm > 0) {
523:     for (i=0;i<lKZ;i++) PetscArraycpy(&data->XKZ[i*data->ldXKZ+i],&data->XKZ[(i+rm)*data->ldXKZ+i+rm],lKZ);
524:   }
525:   lKZ = PetscMin(0,lKZ+V_new);
526:   BVSetActiveColumns(data->KZ,lKZ,lKZ+n);
527:   BVSetActiveColumns(data->U,lKZ,lKZ+n);

529:   /* Compute X, KZ and KR */
530:   BVGetColumn(data->U,lKZ,u);
531:   if (n>1) BVGetColumn(data->U,lKZ+1,&u[1]);
532:   BVGetColumn(data->KZ,lKZ,v);
533:   if (n>1) BVGetColumn(data->KZ,lKZ+1,&v[1]);
534:   d->improvex_jd_proj_uv(d,i_s,i_e,u,v,kr,theta,thetai,pX,pY,ld);
535:   BVRestoreColumn(data->U,lKZ,u);
536:   if (n>1) BVRestoreColumn(data->U,lKZ+1,&u[1]);
537:   BVRestoreColumn(data->KZ,lKZ,v);
538:   if (n>1) BVRestoreColumn(data->KZ,lKZ+1,&v[1]);

540:   /* XKZ <- U'*KZ */
541:   MatCreateSeqDense(PETSC_COMM_SELF,lKZ+n,lKZ+n,NULL,&M);
542:   BVMatProject(data->KZ,NULL,data->U,M);
543:   MatDenseGetArrayRead(M,&array);
544:   for (i=lKZ;i<lKZ+n;i++) { /* upper part */
545:     PetscArraycpy(&data->XKZ[data->ldXKZ*i],&array[i*(lKZ+n)],lKZ);
546:   }
547:   for (i=0;i<lKZ+n;i++) { /* lower part */
548:     PetscArraycpy(&data->XKZ[data->ldXKZ*i+lKZ],&array[i*(lKZ+n)+lKZ],n);
549:   }
550:   MatDenseRestoreArrayRead(M,&array);
551:   MatDestroy(&M);

553:   /* iXKZ <- inv(XKZ) */
554:   size_KZ = lKZ+n;
555:   PetscBLASIntCast(lKZ+n,&s);
556:   data->ldiXKZ = data->size_iXKZ = size_KZ;
557:   for (i=0;i<size_KZ;i++) PetscArraycpy(&data->iXKZ[data->ldiXKZ*i],&data->XKZ[data->ldXKZ*i],size_KZ);
558:   PetscBLASIntCast(data->ldiXKZ,&ldXKZ);
559:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
560:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&s,&s,data->iXKZ,&ldXKZ,data->iXKZPivots,&info));
561:   PetscFPTrapPop();
562:   SlepcCheckLapackInfo("getrf",info);
563:   PetscFunctionReturn(0);
564: }

566: static PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
567: {
568:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
569:   PetscInt       i,j,n,maxits,maxits0,lits,s,ld,k,max_size_D,lV,kV;
570:   PetscScalar    *pX,*pY;
571:   PetscReal      tol,tol0;
572:   Vec            *kr,kr_comp,D_comp,D[2],kr0[2];
573:   PetscBool      odd_situation = PETSC_FALSE;

575:   BVGetActiveColumns(d->eps->V,&lV,&kV);
576:   max_size_D = d->eps->ncv-kV;
577:   /* Quick exit */
578:   if ((max_size_D == 0) || r_e-r_s <= 0) {
579:    *size_D = 0;
580:     PetscFunctionReturn(0);
581:   }

583:   n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);
584:   PetscAssert(n>0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"n == 0");
585:   PetscAssert(data->size_X>=r_e-r_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"size_X < r_e-r_s");

587:   DSGetLeadingDimension(d->eps->ds,&ld);

589:   /* Restart lastTol if a new pair converged */
590:   if (data->dynamic && data->size_cX < lV)
591:     data->lastTol = 0.5;

593:   for (i=0,s=0;i<n;i+=s) {
594:     /* If the selected eigenvalue is complex, but the arithmetic is real... */
595: #if !defined(PETSC_USE_COMPLEX)
596:     if (d->eigi[r_s+i] != 0.0) {
597:       if (i+2 <= max_size_D) s=2;
598:       else break;
599:     } else
600: #endif
601:       s=1;

603:     data->r_s = r_s+i;
604:     data->r_e = r_s+i+s;
605:     SlepcVecPoolGetVecs(d->auxV,s,&kr);

607:     /* Compute theta, maximum iterations and tolerance */
608:     maxits = 0;
609:     tol = 1;
610:     for (j=0;j<s;j++) {
611:       d->improvex_jd_lit(d,r_s+i+j,&data->theta[2*j],&data->thetai[j],&maxits0,&tol0);
612:       maxits += maxits0;
613:       tol *= tol0;
614:     }
615:     maxits/= s;
616:     tol = data->dynamic?data->lastTol:PetscExpReal(PetscLogReal(tol)/s);

618:     /* Compute u, v and kr */
619:     k = r_s+i;
620:     DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);
621:     k = r_s+i;
622:     DSVectors(d->eps->ds,DS_MAT_Y,&k,NULL);
623:     DSGetArray(d->eps->ds,DS_MAT_X,&pX);
624:     DSGetArray(d->eps->ds,DS_MAT_Y,&pY);
625:     dvd_improvex_jd_proj_cuv(d,r_s+i,r_s+i+s,kr,data->theta,data->thetai,pX,pY,ld);
626:     DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);
627:     DSRestoreArray(d->eps->ds,DS_MAT_Y,&pY);

629:     /* Check if the first eigenpairs are converged */
630:     if (i == 0) {
631:       PetscInt oldnpreconv = d->npreconv;
632:       d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv);
633:       if (d->npreconv > oldnpreconv) break;
634:     }

636:     /* Test the odd situation of solving Ax=b with A=I */
637: #if !defined(PETSC_USE_COMPLEX)
638:     odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && data->thetai[0] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
639: #else
640:     odd_situation = (data->ksp && data->theta[0] == 1. && data->theta[1] == 0. && d->B == NULL)? PETSC_TRUE: PETSC_FALSE;
641: #endif
642:     /* If JD */
643:     if (data->ksp && !odd_situation) {
644:       /* kr <- -kr */
645:       for (j=0;j<s;j++) VecScale(kr[j],-1.0);

647:       /* Compose kr and D */
648:       kr0[0] = kr[0];
649:       kr0[1] = (s==2 ? kr[1] : NULL);
650:       VecCreateCompWithVecs(kr0,data->ksp_max_size,data->friends,&kr_comp);
651:       BVGetColumn(d->eps->V,kV+i,&D[0]);
652:       if (s==2) BVGetColumn(d->eps->V,kV+i+1,&D[1]);
653:       else D[1] = NULL;
654:       VecCreateCompWithVecs(D,data->ksp_max_size,data->friends,&D_comp);
655:       VecCompSetSubVecs(data->friends,s,NULL);

657:       /* Solve the correction equation */
658:       KSPSetTolerances(data->ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,maxits);
659:       KSPSolve(data->ksp,kr_comp,D_comp);
660:       KSPGetIterationNumber(data->ksp,&lits);

662:       /* Destroy the composed ks and D */
663:       VecDestroy(&kr_comp);
664:       VecDestroy(&D_comp);
665:       BVRestoreColumn(d->eps->V,kV+i,&D[0]);
666:       if (s==2) BVRestoreColumn(d->eps->V,kV+i+1,&D[1]);

668:     /* If GD */
669:     } else {
670:       BVGetColumn(d->eps->V,kV+i,&D[0]);
671:       if (s==2) BVGetColumn(d->eps->V,kV+i+1,&D[1]);
672:       for (j=0;j<s;j++) d->improvex_precond(d,r_s+i+j,kr[j],D[j]);
673:       dvd_improvex_apply_proj(d,D,s);
674:       BVRestoreColumn(d->eps->V,kV+i,&D[0]);
675:       if (s==2) BVRestoreColumn(d->eps->V,kV+i+1,&D[1]);
676:     }
677:     /* Prevent that short vectors are discarded in the orthogonalization */
678:     if (i == 0 && d->eps->errest[d->nconv+r_s] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s] < PETSC_MAX_REAL) {
679:       for (j=0;j<s;j++) BVScaleColumn(d->eps->V,kV+i+j,1.0/d->eps->errest[d->nconv+r_s]);
680:     }
681:     SlepcVecPoolRestoreVecs(d->auxV,s,&kr);
682:   }
683:   *size_D = i;
684:   if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);
685:   PetscFunctionReturn(0);
686: }

688: PetscErrorCode dvd_improvex_jd(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs,PetscBool dynamic)
689: {
690:   dvdImprovex_jd *data;
691:   PetscBool      useGD;
692:   PC             pc;
693:   PetscInt       size_P;

695:   /* Setting configuration constrains */
696:   PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&useGD);

698:   /* If the arithmetic is real and the problem is not Hermitian, then
699:      the block size is incremented in one */
700: #if !defined(PETSC_USE_COMPLEX)
701:   if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) {
702:     max_bs++;
703:     b->max_size_P = PetscMax(b->max_size_P,2);
704:   } else
705: #endif
706:   {
707:     b->max_size_P = PetscMax(b->max_size_P,1);
708:   }
709:   b->max_size_X = PetscMax(b->max_size_X,max_bs);
710:   size_P = b->max_size_P;

712:   /* Setup the preconditioner */
713:   if (ksp) {
714:     KSPGetPC(ksp,&pc);
715:     dvd_static_precond_PC(d,b,pc);
716:   } else dvd_static_precond_PC(d,b,0);

718:   /* Setup the step */
719:   if (b->state >= DVD_STATE_CONF) {
720:     PetscNewLog(d->eps,&data);
721:     data->dynamic = dynamic;
722:     PetscMalloc1(size_P*size_P,&data->XKZ);
723:     PetscMalloc1(size_P*size_P,&data->iXKZ);
724:     PetscMalloc1(size_P,&data->iXKZPivots);
725:     data->ldXKZ = size_P;
726:     data->size_X = b->max_size_X;
727:     d->improveX_data = data;
728:     data->ksp = useGD? NULL: ksp;
729:     data->d = d;
730:     d->improveX = dvd_improvex_jd_gen;
731: #if !defined(PETSC_USE_COMPLEX)
732:     if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) data->ksp_max_size = 2;
733:     else
734: #endif
735:       data->ksp_max_size = 1;
736:     /* Create various vector basis */
737:     BVDuplicateResize(d->eps->V,size_P,&data->KZ);
738:     BVSetMatrix(data->KZ,NULL,PETSC_FALSE);
739:     BVDuplicate(data->KZ,&data->U);

741:     EPSDavidsonFLAdd(&d->startList,dvd_improvex_jd_start);
742:     EPSDavidsonFLAdd(&d->endList,dvd_improvex_jd_end);
743:     EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_jd_d);
744:   }
745:   PetscFunctionReturn(0);
746: }

748: #if !defined(PETSC_USE_COMPLEX)
749: static inline PetscErrorCode dvd_complex_rayleigh_quotient(Vec ur,Vec ui,Vec Axr,Vec Axi,Vec Bxr,Vec Bxi,PetscScalar *eigr,PetscScalar *eigi)
750: {
751:   PetscScalar    rAr,iAr,rAi,iAi,rBr,iBr,rBi,iBi,b0,b2,b4,b6,b7;

753:   /* eigr = [(rAr+iAi)*(rBr+iBi) + (rAi-iAr)*(rBi-iBr)]/k
754:      eigi = [(rAi-iAr)*(rBr+iBi) - (rAr+iAi)*(rBi-iBr)]/k
755:      k    =  (rBr+iBi)*(rBr+iBi) + (rBi-iBr)*(rBi-iBr)    */
756:   VecDotBegin(Axr,ur,&rAr); /* r*A*r */
757:   VecDotBegin(Axr,ui,&iAr); /* i*A*r */
758:   VecDotBegin(Axi,ur,&rAi); /* r*A*i */
759:   VecDotBegin(Axi,ui,&iAi); /* i*A*i */
760:   VecDotBegin(Bxr,ur,&rBr); /* r*B*r */
761:   VecDotBegin(Bxr,ui,&iBr); /* i*B*r */
762:   VecDotBegin(Bxi,ur,&rBi); /* r*B*i */
763:   VecDotBegin(Bxi,ui,&iBi); /* i*B*i */
764:   VecDotEnd(Axr,ur,&rAr); /* r*A*r */
765:   VecDotEnd(Axr,ui,&iAr); /* i*A*r */
766:   VecDotEnd(Axi,ur,&rAi); /* r*A*i */
767:   VecDotEnd(Axi,ui,&iAi); /* i*A*i */
768:   VecDotEnd(Bxr,ur,&rBr); /* r*B*r */
769:   VecDotEnd(Bxr,ui,&iBr); /* i*B*r */
770:   VecDotEnd(Bxi,ur,&rBi); /* r*B*i */
771:   VecDotEnd(Bxi,ui,&iBi); /* i*B*i */
772:   b0 = rAr+iAi; /* rAr+iAi */
773:   b2 = rAi-iAr; /* rAi-iAr */
774:   b4 = rBr+iBi; /* rBr+iBi */
775:   b6 = rBi-iBr; /* rBi-iBr */
776:   b7 = b4*b4 + b6*b6; /* k */
777:   *eigr = (b0*b4 + b2*b6) / b7; /* eig_r */
778:   *eigi = (b2*b4 - b0*b6) / b7; /* eig_i */
779:   PetscFunctionReturn(0);
780: }
781: #endif

783: static inline PetscErrorCode dvd_compute_n_rr(PetscInt i_s,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,Vec *u,Vec *Ax,Vec *Bx)
784: {
785:   PetscInt       i;
786:   PetscScalar    b0,b1;

788:   for (i=0; i<n; i++) {
789: #if !defined(PETSC_USE_COMPLEX)
790:     if (eigi[i_s+i] != 0.0) {
791:       PetscScalar eigr0=0.0,eigi0=0.0;
792:       dvd_complex_rayleigh_quotient(u[i],u[i+1],Ax[i],Ax[i+1],Bx[i],Bx[i+1],&eigr0,&eigi0);
793:       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) 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);
794:       i++;
795:     } else
796: #endif
797:     {
798:       VecDotBegin(Ax[i],u[i],&b0);
799:       VecDotBegin(Bx[i],u[i],&b1);
800:       VecDotEnd(Ax[i],u[i],&b0);
801:       VecDotEnd(Bx[i],u[i],&b1);
802:       b0 = b0/b1;
803:       if (PetscAbsScalar(eigr[i_s+i]-b0)/PetscAbsScalar(eigr[i_s+i]) > 1e-10) 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));
804:     }
805:   }
806:   PetscFunctionReturn(0);
807: }

809: /*
810:   Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,
811:   kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]
812:   where
813:   pX,pY, the right and left eigenvectors of the projected system
814:   ld, the leading dimension of pX and pY
815: */
816: 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)
817: {
818:   PetscInt       n = i_e-i_s,i;
819:   PetscScalar    *b;
820:   Vec            *Ax,*Bx,*r;
821:   Mat            M;
822:   BV             X;

824:   BVDuplicateResize(d->eps->V,4,&X);
825:   MatCreateSeqDense(PETSC_COMM_SELF,4,4,NULL,&M);
826:   /* u <- X(i) */
827:   dvd_improvex_compute_X(d,i_s,i_e,u,pX,ld);

829:   /* v <- theta[0]A*u + theta[1]*B*u */

831:   /* Bx <- B*X(i) */
832:   Bx = kr;
833:   if (d->BX) {
834:     for (i=i_s; i<i_e; ++i) BVMultVec(d->BX,1.0,0.0,Bx[i-i_s],&pX[ld*i]);
835:   } else {
836:     for (i=0;i<n;i++) {
837:       if (d->B) MatMult(d->B, u[i], Bx[i]);
838:       else VecCopy(u[i], Bx[i]);
839:     }
840:   }

842:   /* Ax <- A*X(i) */
843:   SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&r);
844:   Ax = r;
845:   for (i=i_s; i<i_e; ++i) BVMultVec(d->AX,1.0,0.0,Ax[i-i_s],&pX[ld*i]);

847:   /* v <- Y(i) */
848:   for (i=i_s; i<i_e; ++i) BVMultVec(d->W?d->W:d->eps->V,1.0,0.0,v[i-i_s],&pY[ld*i]);

850:   /* Recompute the eigenvalue */
851:   dvd_compute_n_rr(i_s,n,d->eigr,d->eigi,v,Ax,Bx);

853:   for (i=0;i<n;i++) {
854: #if !defined(PETSC_USE_COMPLEX)
855:     if (d->eigi[i_s+i] != 0.0) {
856:       /* [r_i r_i+1 kr_i kr_i+1]*= [ theta_2i'    0            1        0
857:                                        0         theta_2i'     0        1
858:                                      theta_2i+1 -thetai_i   -eigr_i -eigi_i
859:                                      thetai_i    theta_2i+1  eigi_i -eigr_i ] */
860:       MatDenseGetArrayWrite(M,&b);
861:       b[0] = b[5] = PetscConj(theta[2*i]);
862:       b[2] = b[7] = -theta[2*i+1];
863:       b[6] = -(b[3] = thetai[i]);
864:       b[1] = b[4] = 0.0;
865:       b[8] = b[13] = 1.0/d->nX[i_s+i];
866:       b[10] = b[15] = -d->eigr[i_s+i]/d->nX[i_s+i];
867:       b[14] = -(b[11] = d->eigi[i_s+i]/d->nX[i_s+i]);
868:       b[9] = b[12] = 0.0;
869:       MatDenseRestoreArrayWrite(M,&b);
870:       BVInsertVec(X,0,Ax[i]);
871:       BVInsertVec(X,1,Ax[i+1]);
872:       BVInsertVec(X,2,Bx[i]);
873:       BVInsertVec(X,3,Bx[i+1]);
874:       BVSetActiveColumns(X,0,4);
875:       BVMultInPlace(X,M,0,4);
876:       BVCopyVec(X,0,Ax[i]);
877:       BVCopyVec(X,1,Ax[i+1]);
878:       BVCopyVec(X,2,Bx[i]);
879:       BVCopyVec(X,3,Bx[i+1]);
880:       i++;
881:     } else
882: #endif
883:     {
884:       /* [Ax_i Bx_i]*= [ theta_2i'    1/nX_i
885:                         theta_2i+1  -eig_i/nX_i ] */
886:       MatDenseGetArrayWrite(M,&b);
887:       b[0] = PetscConj(theta[i*2]);
888:       b[1] = theta[i*2+1];
889:       b[4] = 1.0/d->nX[i_s+i];
890:       b[5] = -d->eigr[i_s+i]/d->nX[i_s+i];
891:       MatDenseRestoreArrayWrite(M,&b);
892:       BVInsertVec(X,0,Ax[i]);
893:       BVInsertVec(X,1,Bx[i]);
894:       BVSetActiveColumns(X,0,2);
895:       BVMultInPlace(X,M,0,2);
896:       BVCopyVec(X,0,Ax[i]);
897:       BVCopyVec(X,1,Bx[i]);
898:     }
899:   }
900:   for (i=0; i<n; i++) d->nX[i_s+i] = 1.0;

902:   /* v <- K^{-1} r = K^{-1}(theta_2i'*Ax + theta_2i+1*Bx) */
903:   for (i=0;i<n;i++) d->improvex_precond(d,i_s+i,r[i],v[i]);
904:   SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&r);

906:   /* kr <- P*(Ax - eig_i*Bx) */
907:   d->calcpairs_proj_res(d,i_s,i_e,kr);
908:   BVDestroy(&X);
909:   MatDestroy(&M);
910:   PetscFunctionReturn(0);
911: }

913: static PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d,PetscInt i,PetscScalar* theta,PetscScalar* thetai,PetscInt *maxits,PetscReal *tol)
914: {
915:   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
916:   PetscReal      a;

918:   a = SlepcAbsEigenvalue(d->eigr[i],d->eigi[i]);

920:   if (d->nR[i] < data->fix*a) {
921:     theta[0] = d->eigr[i];
922:     theta[1] = 1.0;
923: #if !defined(PETSC_USE_COMPLEX)
924:     *thetai = d->eigi[i];
925: #endif
926:   } else {
927:     theta[0] = d->target[0];
928:     theta[1] = d->target[1];
929: #if !defined(PETSC_USE_COMPLEX)
930:     *thetai = 0.0;
931: #endif
932: }

934: #if defined(PETSC_USE_COMPLEX)
935:   if (thetai) *thetai = 0.0;
936: #endif
937:   *maxits = data->maxits;
938:   *tol = data->tol;
939:   PetscFunctionReturn(0);
940: }

942: PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard *d,dvdBlackboard *b,PetscInt maxits,PetscReal tol,PetscReal fix)
943: {
944:   dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

946:   /* Setup the step */
947:   if (b->state >= DVD_STATE_CONF) {
948:     data->maxits = maxits;
949:     data->tol = tol;
950:     data->fix = fix;
951:     d->improvex_jd_lit = dvd_improvex_jd_lit_const_0;
952:   }
953:   PetscFunctionReturn(0);
954: }

956: PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d,dvdBlackboard *b)
957: {
958:   /* Setup the step */
959:   if (b->state >= DVD_STATE_CONF) {
960:     d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX;
961:   }
962:   PetscFunctionReturn(0);
963: }

965: PetscErrorCode dvd_improvex_compute_X(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u_,PetscScalar *pX,PetscInt ld)
966: {
967:   PetscInt       n = i_e - i_s,i;
968:   Vec            *u;

970:   if (u_) u = u_;
971:   else if (d->correctXnorm) SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&u);
972:   if (u_ || d->correctXnorm) {
973:     for (i=0;i<n;i++) BVMultVec(d->eps->V,1.0,0.0,u[i],&pX[ld*(i+i_s)]);
974:   }
975:   /* nX(i) <- ||X(i)|| */
976:   if (d->correctXnorm) {
977:     for (i=0;i<n;i++) VecNormBegin(u[i],NORM_2,&d->nX[i_s+i]);
978:     for (i=0;i<n;i++) VecNormEnd(u[i],NORM_2,&d->nX[i_s+i]);
979: #if !defined(PETSC_USE_COMPLEX)
980:     for (i=0;i<n;i++) {
981:       if (d->eigi[i_s+i] != 0.0) {
982:         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]);
983:         i++;
984:       }
985:     }
986: #endif
987:   } else {
988:     for (i=0;i<n;i++) d->nX[i_s+i] = 1.0;
989:   }
990:   if (d->correctXnorm && !u_) SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&u);
991:   PetscFunctionReturn(0);
992: }