Actual source code: dvdimprovex.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }