LCOV - code coverage report
Current view: top level - eps/impls/davidson - dvdimprovex.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 526 578 91.0 %
Date: 2019-10-20 05:41:02 Functions: 23 23 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
       5             : 
       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"
      12             : 
      13             :    Step: improve the eigenvectors X
      14             : */
      15             : 
      16             : #include "davidson.h"
      17             : #include <slepcblaslapack.h>
      18             : 
      19             : /**** JD update step (I - Kfg'/(g'Kf)) K(A - sB) (I - Kfg'/(g'Kf)) t = (I - Kfg'/(g'Kf))r  ****/
      20             : 
      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;
      46             : 
      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       20816 : static PetscErrorCode dvd_improvex_apply_proj(dvdDashboard *d,Vec *V,PetscInt cV)
      53             : {
      54             : #if defined(PETSC_MISSING_LAPACK_GETRS)
      55             :   PetscFunctionBegin;
      56             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable");
      57             : #else
      58             :   PetscErrorCode ierr;
      59       20816 :   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
      66             : 
      67       20816 :   PetscFunctionBegin;
      68       20816 :   if (cV > 2) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
      69             : 
      70             :   /* h <- X'*V */
      71       20816 :   ierr = PetscMalloc1(data->size_iXKZ*cV,&h);CHKERRQ(ierr);
      72       20816 :   ldh = data->size_iXKZ;
      73       20816 :   ierr = BVGetActiveColumns(data->U,&l,&k);CHKERRQ(ierr);
      74       20816 :   if (ldh!=k) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
      75       20816 :   ierr = BVSetActiveColumns(data->U,0,k);CHKERRQ(ierr);
      76       20873 :   for (i=0;i<cV;i++) {
      77       20873 :     ierr = BVDotVec(data->U,V[i],&h[ldh*i]);CHKERRQ(ierr);
      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       20816 :   ierr = BVSetActiveColumns(data->U,l,k);CHKERRQ(ierr);
      83             : 
      84             :   /* h <- iXKZ\h */
      85       20816 :   ierr = PetscBLASIntCast(cV,&cV_);CHKERRQ(ierr);
      86       20816 :   ierr = PetscBLASIntCast(data->size_iXKZ,&n);CHKERRQ(ierr);
      87       20816 :   ierr = PetscBLASIntCast(data->ldiXKZ,&ld);CHKERRQ(ierr);
      88       20816 :   PetscValidScalarPointer(data->iXKZ,0);
      89       20816 :   ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
      90       20816 :   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
      91       20816 :   ierr = PetscFPTrapPop();CHKERRQ(ierr);
      92       20816 :   SlepcCheckLapackInfo("getrs",info);
      93             : 
      94             :   /* V <- V - KZ*h */
      95       20816 :   ierr = BVSetActiveColumns(data->KZ,0,k);CHKERRQ(ierr);
      96       20873 :   for (i=0;i<cV;i++) {
      97       20873 :     ierr = BVMultVec(data->KZ,-1.0,1.0,V[i],&h[ldh*i]);CHKERRQ(ierr);
      98             :   }
      99       20816 :   ierr = BVSetActiveColumns(data->KZ,l,k);CHKERRQ(ierr);
     100       20816 :   ierr = PetscFree(h);CHKERRQ(ierr);
     101       20816 :   PetscFunctionReturn(0);
     102             : #endif
     103             : }
     104             : 
     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          96 : static PetscErrorCode dvd_improvex_applytrans_proj(dvdDashboard *d,Vec *V,PetscInt cV)
     111             : {
     112             : #if defined(PETSC_MISSING_LAPACK_GETRS)
     113             :   PetscFunctionBegin;
     114             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable");
     115             :   PetscFunctionReturn(0);
     116             : #else
     117             :   PetscErrorCode ierr;
     118          96 :   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
     125             : 
     126          96 :   PetscFunctionBegin;
     127          96 :   if (cV > 2) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     128             : 
     129             :   /* h <- KZ'*V */
     130          96 :   ierr = PetscMalloc1(data->size_iXKZ*cV,&h);CHKERRQ(ierr);
     131          96 :   ldh = data->size_iXKZ;
     132          96 :   ierr = BVGetActiveColumns(data->U,&l,&k);CHKERRQ(ierr);
     133          96 :   if (ldh!=k) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     134          96 :   ierr = BVSetActiveColumns(data->KZ,0,k);CHKERRQ(ierr);
     135          96 :   for (i=0;i<cV;i++) {
     136          96 :     ierr = BVDotVec(data->KZ,V[i],&h[ldh*i]);CHKERRQ(ierr);
     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          96 :   ierr = BVSetActiveColumns(data->KZ,l,k);CHKERRQ(ierr);
     142             : 
     143             :   /* h <- iXKZ\h */
     144          96 :   ierr = PetscBLASIntCast(cV,&cV_);CHKERRQ(ierr);
     145          96 :   ierr = PetscBLASIntCast(data->size_iXKZ,&n);CHKERRQ(ierr);
     146          96 :   ierr = PetscBLASIntCast(data->ldiXKZ,&ld);CHKERRQ(ierr);
     147          96 :   PetscValidScalarPointer(data->iXKZ,0);
     148          96 :   ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
     149          96 :   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&cV_,data->iXKZ,&ld,data->iXKZPivots,h,&n,&info));
     150          96 :   ierr = PetscFPTrapPop();CHKERRQ(ierr);
     151          96 :   SlepcCheckLapackInfo("getrs",info);
     152             : 
     153             :   /* V <- V - U*h */
     154          96 :   ierr = BVSetActiveColumns(data->U,0,k);CHKERRQ(ierr);
     155          96 :   for (i=0;i<cV;i++) {
     156          96 :     ierr = BVMultVec(data->U,-1.0,1.0,V[i],&h[ldh*i]);CHKERRQ(ierr);
     157             :   }
     158          96 :   ierr = BVSetActiveColumns(data->U,l,k);CHKERRQ(ierr);
     159          96 :   ierr = PetscFree(h);CHKERRQ(ierr);
     160          96 :   PetscFunctionReturn(0);
     161             : #endif
     162             : }
     163             : 
     164          73 : static PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)
     165             : {
     166             :   PetscErrorCode ierr;
     167          73 :   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
     168             : 
     169          73 :   PetscFunctionBegin;
     170          73 :   ierr = VecDestroy(&data->friends);CHKERRQ(ierr);
     171             : 
     172             :   /* Restore the pc of ksp */
     173          73 :   if (data->old_pc) {
     174          19 :     ierr = KSPSetPC(data->ksp, data->old_pc);CHKERRQ(ierr);
     175          19 :     ierr = PCDestroy(&data->old_pc);CHKERRQ(ierr);
     176             :   }
     177          73 :   PetscFunctionReturn(0);
     178             : }
     179             : 
     180          73 : static PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)
     181             : {
     182             :   PetscErrorCode ierr;
     183          73 :   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
     184             : 
     185          73 :   PetscFunctionBegin;
     186             :   /* Free local data and objects */
     187          73 :   ierr = PetscFree(data->XKZ);CHKERRQ(ierr);
     188          73 :   ierr = PetscFree(data->iXKZ);CHKERRQ(ierr);
     189          73 :   ierr = PetscFree(data->iXKZPivots);CHKERRQ(ierr);
     190          73 :   ierr = BVDestroy(&data->KZ);CHKERRQ(ierr);
     191          73 :   ierr = BVDestroy(&data->U);CHKERRQ(ierr);
     192          73 :   ierr = PetscFree(data);CHKERRQ(ierr);
     193          73 :   PetscFunctionReturn(0);
     194             : }
     195             : 
     196             : /*
     197             :    y <- theta[1]A*x - theta[0]*B*x
     198             :    auxV, two auxiliary vectors
     199             :  */
     200       29089 : PETSC_STATIC_INLINE PetscErrorCode dvd_aux_matmult(dvdImprovex_jd *data,const Vec *x,const Vec *y)
     201             : {
     202             :   PetscErrorCode ierr;
     203             :   PetscInt       n,i;
     204             :   const Vec      *Bx;
     205             :   Vec            *auxV;
     206             : 
     207       29089 :   PetscFunctionBegin;
     208       29089 :   n = data->r_e - data->r_s;
     209       58178 :   for (i=0;i<n;i++) {
     210       29089 :     ierr = MatMult(data->d->A,x[i],y[i]);CHKERRQ(ierr);
     211             :   }
     212             : 
     213       29089 :   ierr = SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);CHKERRQ(ierr);
     214       29089 :   for (i=0;i<n;i++) {
     215             : #if !defined(PETSC_USE_COMPLEX)
     216       29089 :     if (data->d->eigi[data->r_s+i] != 0.0) {
     217           0 :       if (data->d->B) {
     218           0 :         ierr = MatMult(data->d->B,x[i],auxV[0]);CHKERRQ(ierr);
     219           0 :         ierr = MatMult(data->d->B,x[i+1],auxV[1]);CHKERRQ(ierr);
     220           0 :         Bx = auxV;
     221           0 :       } else Bx = &x[i];
     222             : 
     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           0 :       ierr = VecAXPBYPCZ(y[i],-data->theta[2*i],data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);CHKERRQ(ierr);
     226           0 :       ierr = VecAXPBYPCZ(y[i+1],-data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);CHKERRQ(ierr);
     227           0 :       i++;
     228             :     } else
     229             : #endif
     230             :     {
     231       29089 :       if (data->d->B) {
     232        7864 :         ierr = MatMult(data->d->B,x[i],auxV[0]);CHKERRQ(ierr);
     233        7864 :         Bx = auxV;
     234       21225 :       } else Bx = &x[i];
     235       29089 :       ierr = VecAXPBY(y[i],-data->theta[i*2],data->theta[i*2+1],Bx[0]);CHKERRQ(ierr);
     236             :     }
     237             :   }
     238       29089 :   ierr = SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);CHKERRQ(ierr);
     239       29089 :   PetscFunctionReturn(0);
     240             : }
     241             : 
     242             : /*
     243             :    y <- theta[1]'*A'*x - theta[0]'*B'*x
     244             :  */
     245          83 : PETSC_STATIC_INLINE PetscErrorCode dvd_aux_matmulttrans(dvdImprovex_jd *data,const Vec *x,const Vec *y)
     246             : {
     247             :   PetscErrorCode ierr;
     248             :   PetscInt       n,i;
     249             :   const Vec      *Bx;
     250             :   Vec            *auxV;
     251             : 
     252          83 :   PetscFunctionBegin;
     253          83 :   n = data->r_e - data->r_s;
     254         166 :   for (i=0;i<n;i++) {
     255          83 :     ierr = MatMultTranspose(data->d->A,x[i],y[i]);CHKERRQ(ierr);
     256             :   }
     257             : 
     258          83 :   ierr = SlepcVecPoolGetVecs(data->d->auxV,2,&auxV);CHKERRQ(ierr);
     259          83 :   for (i=0;i<n;i++) {
     260             : #if !defined(PETSC_USE_COMPLEX)
     261          83 :     if (data->d->eigi[data->r_s+i] != 0.0) {
     262           0 :       if (data->d->B) {
     263           0 :         ierr = MatMultTranspose(data->d->B,x[i],auxV[0]);CHKERRQ(ierr);
     264           0 :         ierr = MatMultTranspose(data->d->B,x[i+1],auxV[1]);CHKERRQ(ierr);
     265           0 :         Bx = auxV;
     266           0 :       } else Bx = &x[i];
     267             : 
     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           0 :       ierr = VecAXPBYPCZ(y[i],-data->theta[2*i],-data->thetai[i],data->theta[2*i+1],Bx[0],Bx[1]);CHKERRQ(ierr);
     271           0 :       ierr = VecAXPBYPCZ(y[i+1],data->thetai[i],-data->theta[2*i],data->theta[2*i+1],Bx[0],Bx[1]);CHKERRQ(ierr);
     272           0 :       i++;
     273             :     } else
     274             : #endif
     275             :     {
     276          83 :       if (data->d->B) {
     277           0 :         ierr = MatMultTranspose(data->d->B,x[i],auxV[0]);CHKERRQ(ierr);
     278           0 :         Bx = auxV;
     279          83 :       } else Bx = &x[i];
     280          83 :       ierr = VecAXPBY(y[i],PetscConj(-data->theta[i*2]),PetscConj(data->theta[i*2+1]),Bx[0]);CHKERRQ(ierr);
     281             :     }
     282             :   }
     283          83 :   ierr = SlepcVecPoolRestoreVecs(data->d->auxV,2,&auxV);CHKERRQ(ierr);
     284          83 :   PetscFunctionReturn(0);
     285             : }
     286             : 
     287       15596 : static PetscErrorCode PCApplyBA_dvd(PC pc,PCSide side,Vec in,Vec out,Vec w)
     288             : {
     289             :   PetscErrorCode ierr;
     290             :   dvdImprovex_jd *data;
     291             :   PetscInt       n,i;
     292             :   const Vec      *inx,*outx,*wx;
     293             :   Vec            *auxV;
     294             :   Mat            A;
     295             : 
     296       15596 :   PetscFunctionBegin;
     297       15596 :   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr);
     298       15596 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     299       15596 :   ierr = VecCompGetSubVecs(in,NULL,&inx);CHKERRQ(ierr);
     300       15596 :   ierr = VecCompGetSubVecs(out,NULL,&outx);CHKERRQ(ierr);
     301       15596 :   ierr = VecCompGetSubVecs(w,NULL,&wx);CHKERRQ(ierr);
     302       15596 :   n = data->r_e - data->r_s;
     303       15596 :   ierr = SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     304       15596 :   switch (side) {
     305             :   case PC_LEFT:
     306             :     /* aux <- theta[1]A*in - theta[0]*B*in */
     307       15596 :     ierr = dvd_aux_matmult(data,inx,auxV);CHKERRQ(ierr);
     308             : 
     309             :     /* out <- K * aux */
     310       15596 :     for (i=0;i<n;i++) {
     311       15596 :       ierr = data->d->improvex_precond(data->d,data->r_s+i,auxV[i],outx[i]);CHKERRQ(ierr);
     312             :     }
     313             :     break;
     314             :   case PC_RIGHT:
     315             :     /* aux <- K * in */
     316           0 :     for (i=0;i<n;i++) {
     317           0 :       ierr = data->d->improvex_precond(data->d,data->r_s+i,inx[i],auxV[i]);CHKERRQ(ierr);
     318             :     }
     319             : 
     320             :     /* out <- theta[1]A*auxV - theta[0]*B*auxV */
     321           0 :     ierr = dvd_aux_matmult(data,auxV,outx);CHKERRQ(ierr);
     322             :     break;
     323             :   case PC_SYMMETRIC:
     324             :     /* aux <- K^{1/2} * in */
     325           0 :     for (i=0;i<n;i++) {
     326           0 :       ierr = PCApplySymmetricRight(data->old_pc,inx[i],auxV[i]);CHKERRQ(ierr);
     327             :     }
     328             : 
     329             :     /* wx <- theta[1]A*auxV - theta[0]*B*auxV */
     330           0 :     ierr = dvd_aux_matmult(data,auxV,wx);CHKERRQ(ierr);
     331             : 
     332             :     /* aux <- K^{1/2} * in */
     333           0 :     for (i=0;i<n;i++) {
     334           0 :       ierr = PCApplySymmetricLeft(data->old_pc,wx[i],outx[i]);CHKERRQ(ierr);
     335             :     }
     336             :     break;
     337             :   default:
     338           0 :     SETERRQ(PETSC_COMM_SELF,1,"Unsupported KSP side");
     339             :   }
     340             :   /* out <- out - v*(u'*out) */
     341       15596 :   ierr = dvd_improvex_apply_proj(data->d,(Vec*)outx,n);CHKERRQ(ierr);
     342       15596 :   ierr = SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     343       15596 :   PetscFunctionReturn(0);
     344             : }
     345             : 
     346         372 : static PetscErrorCode PCApply_dvd(PC pc,Vec in,Vec out)
     347             : {
     348             :   PetscErrorCode ierr;
     349             :   dvdImprovex_jd *data;
     350             :   PetscInt       n,i;
     351             :   const Vec      *inx, *outx;
     352             :   Mat            A;
     353             : 
     354         372 :   PetscFunctionBegin;
     355         372 :   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr);
     356         372 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     357         372 :   ierr = VecCompGetSubVecs(in,NULL,&inx);CHKERRQ(ierr);
     358         372 :   ierr = VecCompGetSubVecs(out,NULL,&outx);CHKERRQ(ierr);
     359         372 :   n = data->r_e - data->r_s;
     360             :   /* out <- K * in */
     361         744 :   for (i=0;i<n;i++) {
     362         372 :     ierr = data->d->improvex_precond(data->d,data->r_s+i,inx[i],outx[i]);CHKERRQ(ierr);
     363             :   }
     364             :   /* out <- out - v*(u'*out) */
     365         372 :   ierr = dvd_improvex_apply_proj(data->d,(Vec*)outx,n);CHKERRQ(ierr);
     366         372 :   PetscFunctionReturn(0);
     367             : }
     368             : 
     369          96 : static PetscErrorCode PCApplyTranspose_dvd(PC pc,Vec in,Vec out)
     370             : {
     371             :   PetscErrorCode ierr;
     372             :   dvdImprovex_jd *data;
     373             :   PetscInt       n,i;
     374             :   const Vec      *inx, *outx;
     375             :   Vec            *auxV;
     376             :   Mat            A;
     377             : 
     378          96 :   PetscFunctionBegin;
     379          96 :   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr);
     380          96 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     381          96 :   ierr = VecCompGetSubVecs(in,NULL,&inx);CHKERRQ(ierr);
     382          96 :   ierr = VecCompGetSubVecs(out,NULL,&outx);CHKERRQ(ierr);
     383          96 :   n = data->r_e - data->r_s;
     384          96 :   ierr = SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     385             :   /* auxV <- in */
     386          96 :   for (i=0;i<n;i++) {
     387          96 :     ierr = VecCopy(inx[i],auxV[i]);CHKERRQ(ierr);
     388             :   }
     389             :   /* auxV <- auxV - u*(v'*auxV) */
     390          96 :   ierr = dvd_improvex_applytrans_proj(data->d,auxV,n);CHKERRQ(ierr);
     391             :   /* out <- K' * aux */
     392          96 :   for (i=0;i<n;i++) {
     393          96 :     ierr = PCApplyTranspose(data->old_pc,auxV[i],outx[i]);CHKERRQ(ierr);
     394             :   }
     395          96 :   ierr = SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     396          96 :   PetscFunctionReturn(0);
     397             : }
     398             : 
     399       13493 : static PetscErrorCode MatMult_dvd_jd(Mat A,Vec in,Vec out)
     400             : {
     401             :   PetscErrorCode ierr;
     402             :   dvdImprovex_jd *data;
     403             :   PetscInt       n;
     404             :   const Vec      *inx, *outx;
     405             :   PCSide         side;
     406             : 
     407       13493 :   PetscFunctionBegin;
     408       13493 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     409       13493 :   ierr = VecCompGetSubVecs(in,NULL,&inx);CHKERRQ(ierr);
     410       13493 :   ierr = VecCompGetSubVecs(out,NULL,&outx);CHKERRQ(ierr);
     411       13493 :   n = data->r_e - data->r_s;
     412             :   /* out <- theta[1]A*in - theta[0]*B*in */
     413       13493 :   ierr = dvd_aux_matmult(data,inx,outx);CHKERRQ(ierr);
     414       13493 :   ierr = KSPGetPCSide(data->ksp,&side);CHKERRQ(ierr);
     415       13493 :   if (side == PC_RIGHT) {
     416             :     /* out <- out - v*(u'*out) */
     417           0 :     ierr = dvd_improvex_apply_proj(data->d,(Vec*)outx,n);CHKERRQ(ierr);
     418             :   }
     419       13493 :   PetscFunctionReturn(0);
     420             : }
     421             : 
     422          83 : static PetscErrorCode MatMultTranspose_dvd_jd(Mat A,Vec in,Vec out)
     423             : {
     424             :   PetscErrorCode ierr;
     425             :   dvdImprovex_jd *data;
     426             :   PetscInt       n,i;
     427             :   const Vec      *inx,*outx,*r;
     428             :   Vec            *auxV;
     429             :   PCSide         side;
     430             : 
     431          83 :   PetscFunctionBegin;
     432          83 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     433          83 :   ierr = VecCompGetSubVecs(in,NULL,&inx);CHKERRQ(ierr);
     434          83 :   ierr = VecCompGetSubVecs(out,NULL,&outx);CHKERRQ(ierr);
     435          83 :   n = data->r_e - data->r_s;
     436          83 :   ierr = KSPGetPCSide(data->ksp,&side);CHKERRQ(ierr);
     437          83 :   if (side == PC_RIGHT) {
     438             :     /* auxV <- in */
     439           0 :     ierr = SlepcVecPoolGetVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     440           0 :     for (i=0;i<n;i++) {
     441           0 :       ierr = VecCopy(inx[i],auxV[i]);CHKERRQ(ierr);
     442             :     }
     443             :     /* auxV <- auxV - v*(u'*auxV) */
     444           0 :     ierr = dvd_improvex_applytrans_proj(data->d,auxV,n);CHKERRQ(ierr);
     445           0 :     r = auxV;
     446          83 :   } else r = inx;
     447             :   /* out <- theta[1]A*r - theta[0]*B*r */
     448          83 :   ierr = dvd_aux_matmulttrans(data,r,outx);CHKERRQ(ierr);
     449          83 :   if (side == PC_RIGHT) {
     450           0 :     ierr = SlepcVecPoolRestoreVecs(data->d->auxV,n,&auxV);CHKERRQ(ierr);
     451             :   }
     452          83 :   PetscFunctionReturn(0);
     453             : }
     454             : 
     455          23 : static PetscErrorCode MatCreateVecs_dvd_jd(Mat A,Vec *right,Vec *left)
     456             : {
     457             :   PetscErrorCode ierr;
     458             :   Vec            *r,*l;
     459             :   dvdImprovex_jd *data;
     460             :   PetscInt       n,i;
     461             : 
     462          23 :   PetscFunctionBegin;
     463          23 :   ierr = MatShellGetContext(A,(void**)&data);CHKERRQ(ierr);
     464          23 :   n = data->ksp_max_size;
     465          23 :   if (right) {
     466          23 :     ierr = PetscMalloc1(n,&r);CHKERRQ(ierr);
     467             :   }
     468          23 :   if (left) {
     469           0 :     ierr = PetscMalloc1(n,&l);CHKERRQ(ierr);
     470             :   }
     471          31 :   for (i=0;i<n;i++) {
     472          31 :     ierr = MatCreateVecs(data->d->A,right?&r[i]:NULL,left?&l[i]:NULL);CHKERRQ(ierr);
     473             :   }
     474          23 :   if (right) {
     475          23 :     ierr = VecCreateCompWithVecs(r,n,data->friends,right);CHKERRQ(ierr);
     476          31 :     for (i=0;i<n;i++) {
     477          31 :       ierr = VecDestroy(&r[i]);CHKERRQ(ierr);
     478             :     }
     479             :   }
     480          23 :   if (left) {
     481           0 :     ierr = VecCreateCompWithVecs(l,n,data->friends,left);CHKERRQ(ierr);
     482           0 :     for (i=0;i<n;i++) {
     483           0 :       ierr = VecDestroy(&l[i]);CHKERRQ(ierr);
     484             :     }
     485             :   }
     486             : 
     487          23 :   if (right) {
     488          23 :     ierr = PetscFree(r);CHKERRQ(ierr);
     489             :   }
     490          23 :   if (left) {
     491           0 :     ierr = PetscFree(l);CHKERRQ(ierr);
     492             :   }
     493          23 :   PetscFunctionReturn(0);
     494             : }
     495             : 
     496          73 : static PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)
     497             : {
     498             :   PetscErrorCode ierr;
     499          73 :   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];
     505             : 
     506          73 :   PetscFunctionBegin;
     507          73 :   data->size_cX = data->old_size_X = 0;
     508          73 :   data->lastTol = data->dynamic?0.5:0.0;
     509             : 
     510             :   /* Setup the ksp */
     511          73 :   if (data->ksp) {
     512             :     /* Create the reference vector */
     513          30 :     ierr = BVGetColumn(d->eps->V,0,&v0[0]);CHKERRQ(ierr);
     514          30 :     v0[1] = v0[0];
     515          30 :     ierr = VecCreateCompWithVecs(v0,data->ksp_max_size,NULL,&data->friends);CHKERRQ(ierr);
     516          30 :     ierr = BVRestoreColumn(d->eps->V,0,&v0[0]);CHKERRQ(ierr);
     517          30 :     ierr = PetscLogObjectParent((PetscObject)d->eps,(PetscObject)data->friends);CHKERRQ(ierr);
     518             : 
     519             :     /* Save the current pc and set a PCNONE */
     520          30 :     ierr = KSPGetPC(data->ksp, &data->old_pc);CHKERRQ(ierr);
     521          30 :     ierr = PetscObjectTypeCompare((PetscObject)data->old_pc,PCNONE,&t);CHKERRQ(ierr);
     522          30 :     data->lastTol = 0.5;
     523          30 :     if (t) data->old_pc = 0;
     524             :     else {
     525          19 :       ierr = PetscObjectReference((PetscObject)data->old_pc);CHKERRQ(ierr);
     526          19 :       ierr = PCCreate(PetscObjectComm((PetscObject)d->eps),&pc);CHKERRQ(ierr);
     527          19 :       ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
     528          19 :       ierr = PCSetOperators(pc,d->A,d->A);CHKERRQ(ierr);
     529          19 :       ierr = PCSetReusePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr);
     530          19 :       ierr = PCShellSetApply(pc,PCApply_dvd);CHKERRQ(ierr);
     531          19 :       ierr = PCShellSetApplyBA(pc,PCApplyBA_dvd);CHKERRQ(ierr);
     532          19 :       ierr = PCShellSetApplyTranspose(pc,PCApplyTranspose_dvd);CHKERRQ(ierr);
     533          19 :       ierr = KSPSetPC(data->ksp,pc);CHKERRQ(ierr);
     534          19 :       ierr = PCDestroy(&pc);CHKERRQ(ierr);
     535             :     }
     536             : 
     537             :     /* Create the (I-v*u')*K*(A-s*B) matrix */
     538          30 :     ierr = MatGetSize(d->A,&rA,&cA);CHKERRQ(ierr);
     539          30 :     ierr = MatGetLocalSize(d->A,&rlA,&clA);CHKERRQ(ierr);
     540          30 :     ierr = 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);CHKERRQ(ierr);
     541          30 :     ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_dvd_jd);CHKERRQ(ierr);
     542          30 :     ierr = MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_dvd_jd);CHKERRQ(ierr);
     543          30 :     ierr = MatShellSetOperation(A,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_dvd_jd);CHKERRQ(ierr);
     544             : 
     545             :     /* Try to avoid KSPReset */
     546          30 :     ierr = KSPGetOperatorsSet(data->ksp,&t,NULL);CHKERRQ(ierr);
     547          30 :     if (t) {
     548             :       Mat      M;
     549             :       PetscInt rM;
     550          30 :       ierr = KSPGetOperators(data->ksp,&M,NULL);CHKERRQ(ierr);
     551          30 :       ierr = MatGetSize(M,&rM,NULL);CHKERRQ(ierr);
     552          30 :       if (rM != rA*data->ksp_max_size) { ierr = KSPReset(data->ksp);CHKERRQ(ierr); }
     553             :     }
     554          30 :     ierr = KSPSetOperators(data->ksp,A,A);CHKERRQ(ierr);
     555          30 :     ierr = KSPSetReusePreconditioner(data->ksp,PETSC_TRUE);CHKERRQ(ierr);
     556          30 :     ierr = KSPSetUp(data->ksp);CHKERRQ(ierr);
     557          30 :     ierr = MatDestroy(&A);CHKERRQ(ierr);
     558             :   } else {
     559          43 :     data->old_pc = 0;
     560          43 :     data->friends = NULL;
     561             :   }
     562          73 :   ierr = BVSetActiveColumns(data->KZ,0,0);CHKERRQ(ierr);
     563          73 :   ierr = BVSetActiveColumns(data->U,0,0);CHKERRQ(ierr);
     564          73 :   PetscFunctionReturn(0);
     565             : }
     566             : 
     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        5562 : 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)
     577             :   PetscFunctionBegin;
     578             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable");
     579             : #else
     580             :   PetscErrorCode ierr;
     581        5562 :   PetscInt       n=i_e-i_s,size_KZ,V_new,rm,i,lv,kv,lKZ,kKZ;
     582        5562 :   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;
     587             : 
     588        5562 :   PetscFunctionBegin;
     589             :   /* Check consistency */
     590        5562 :   ierr = BVGetActiveColumns(d->eps->V,&lv,&kv);CHKERRQ(ierr);
     591        5562 :   V_new = lv - data->size_cX;
     592        5562 :   if (V_new > data->old_size_X) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     593        5562 :   data->old_size_X = n;
     594        5562 :   data->size_cX = lv;
     595             : 
     596             :   /* KZ <- KZ(rm:rm+max_cX-1) */
     597        5562 :   ierr = BVGetActiveColumns(data->KZ,&lKZ,&kKZ);CHKERRQ(ierr);
     598        5562 :   rm = PetscMax(V_new+lKZ,0);
     599        5562 :   if (rm > 0) {
     600           0 :     for (i=0;i<lKZ;i++) {
     601           0 :       ierr = BVCopyColumn(data->KZ,i+rm,i);CHKERRQ(ierr);
     602           0 :       ierr = BVCopyColumn(data->U,i+rm,i);CHKERRQ(ierr);
     603             :     }
     604             :   }
     605             : 
     606             :   /* XKZ <- XKZ(rm:rm+max_cX-1,rm:rm+max_cX-1) */
     607        5562 :   if (rm > 0) {
     608           0 :     for (i=0;i<lKZ;i++) {
     609           0 :       ierr = PetscArraycpy(&data->XKZ[i*data->ldXKZ+i],&data->XKZ[(i+rm)*data->ldXKZ+i+rm],lKZ);CHKERRQ(ierr);
     610             :     }
     611             :   }
     612        5562 :   lKZ = PetscMin(0,lKZ+V_new);
     613        5562 :   ierr = BVSetActiveColumns(data->KZ,lKZ,lKZ+n);CHKERRQ(ierr);
     614        5562 :   ierr = BVSetActiveColumns(data->U,lKZ,lKZ+n);CHKERRQ(ierr);
     615             : 
     616             :   /* Compute X, KZ and KR */
     617        5562 :   ierr = BVGetColumn(data->U,lKZ,u);CHKERRQ(ierr);
     618        5562 :   if (n>1) { ierr = BVGetColumn(data->U,lKZ+1,&u[1]);CHKERRQ(ierr); }
     619        5562 :   ierr = BVGetColumn(data->KZ,lKZ,v);CHKERRQ(ierr);
     620        5562 :   if (n>1) { ierr = BVGetColumn(data->KZ,lKZ+1,&v[1]);CHKERRQ(ierr); }
     621        5562 :   ierr = d->improvex_jd_proj_uv(d,i_s,i_e,u,v,kr,theta,thetai,pX,pY,ld);CHKERRQ(ierr);
     622        5562 :   ierr = BVRestoreColumn(data->U,lKZ,u);CHKERRQ(ierr);
     623        5562 :   if (n>1) { ierr = BVRestoreColumn(data->U,lKZ+1,&u[1]);CHKERRQ(ierr); }
     624        5562 :   ierr = BVRestoreColumn(data->KZ,lKZ,v);CHKERRQ(ierr);
     625        5562 :   if (n>1) { ierr = BVRestoreColumn(data->KZ,lKZ+1,&v[1]);CHKERRQ(ierr); }
     626             : 
     627             :   /* XKZ <- U'*KZ */
     628        5562 :   ierr = MatCreateSeqDense(PETSC_COMM_SELF,lKZ+n,lKZ+n,NULL,&M);CHKERRQ(ierr);
     629        5562 :   ierr = BVMatProject(data->KZ,NULL,data->U,M);CHKERRQ(ierr);
     630        5562 :   ierr = MatDenseGetArray(M,&array);CHKERRQ(ierr);
     631       11183 :   for (i=lKZ;i<lKZ+n;i++) { /* upper part */
     632        5621 :     ierr = PetscArraycpy(&data->XKZ[data->ldXKZ*i],&array[i*(lKZ+n)],lKZ);CHKERRQ(ierr);
     633             :   }
     634        5621 :   for (i=0;i<lKZ+n;i++) { /* lower part */
     635        5621 :     ierr = PetscArraycpy(&data->XKZ[data->ldXKZ*i+lKZ],&array[i*(lKZ+n)+lKZ],n);CHKERRQ(ierr);
     636             :   }
     637        5562 :   ierr = MatDenseRestoreArray(M,&array);CHKERRQ(ierr);
     638        5562 :   ierr = MatDestroy(&M);CHKERRQ(ierr);
     639             : 
     640             :   /* iXKZ <- inv(XKZ) */
     641        5562 :   size_KZ = lKZ+n;
     642        5562 :   ierr = PetscBLASIntCast(lKZ+n,&s);CHKERRQ(ierr);
     643        5562 :   data->ldiXKZ = data->size_iXKZ = size_KZ;
     644       11183 :   for (i=0;i<size_KZ;i++) {
     645        5621 :     ierr = PetscArraycpy(&data->iXKZ[data->ldiXKZ*i],&data->XKZ[data->ldXKZ*i],size_KZ);CHKERRQ(ierr);
     646             :   }
     647        5562 :   ierr = PetscBLASIntCast(data->ldiXKZ,&ldXKZ);CHKERRQ(ierr);
     648        5562 :   ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
     649        5562 :   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&s,&s,data->iXKZ,&ldXKZ,data->iXKZPivots,&info));
     650        5562 :   ierr = PetscFPTrapPop();CHKERRQ(ierr);
     651        5562 :   SlepcCheckLapackInfo("getrf",info);
     652        5562 :   PetscFunctionReturn(0);
     653             : #endif
     654             : }
     655             : 
     656        5025 : static PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
     657             : {
     658        5025 :   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
     659             :   PetscErrorCode ierr;
     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        5025 :   PetscBool      odd_situation = PETSC_FALSE;
     665             : 
     666        5025 :   PetscFunctionBegin;
     667        5025 :   ierr = BVGetActiveColumns(d->eps->V,&lV,&kV);CHKERRQ(ierr);
     668        5025 :   max_size_D = d->eps->ncv-kV;
     669             :   /* Quick exit */
     670        5025 :   if ((max_size_D == 0) || r_e-r_s <= 0) {
     671           0 :    *size_D = 0;
     672           0 :     PetscFunctionReturn(0);
     673             :   }
     674             : 
     675        5025 :   n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);
     676        5025 :   if (n == 0) SETERRQ(PETSC_COMM_SELF,1,"n == 0");
     677        5025 :   if (data->size_X < r_e-r_s) SETERRQ(PETSC_COMM_SELF,1,"size_X < r_e-r_s");
     678             : 
     679        5025 :   ierr = DSGetLeadingDimension(d->eps->ds,&ld);CHKERRQ(ierr);
     680             : 
     681             :   /* Restart lastTol if a new pair converged */
     682        5025 :   if (data->dynamic && data->size_cX < lV)
     683           1 :     data->lastTol = 0.5;
     684             : 
     685        5318 :   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        5563 :     if (d->eigi[r_s+i] != 0.0) {
     689          60 :       if (i+2 <= max_size_D) s=2;
     690             :       else break;
     691             :     } else
     692             : #endif
     693             :       s=1;
     694             : 
     695        5562 :     data->r_s = r_s+i;
     696        5562 :     data->r_e = r_s+i+s;
     697        5562 :     ierr = SlepcVecPoolGetVecs(d->auxV,s,&kr);CHKERRQ(ierr);
     698             : 
     699             :     /* Compute theta, maximum iterations and tolerance */
     700             :     maxits = 0;
     701             :     tol = 1;
     702        5621 :     for (j=0;j<s;j++) {
     703        5621 :       ierr = d->improvex_jd_lit(d,r_s+i+j,&data->theta[2*j],&data->thetai[j],&maxits0,&tol0);CHKERRQ(ierr);
     704        5621 :       maxits += maxits0;
     705        5621 :       tol *= tol0;
     706             :     }
     707        5562 :     maxits/= s;
     708        5562 :     tol = data->dynamic?data->lastTol:PetscExpReal(PetscLogReal(tol)/s);
     709             : 
     710             :     /* Compute u, v and kr */
     711        5562 :     k = r_s+i;
     712        5562 :     ierr = DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);CHKERRQ(ierr);
     713        5562 :     k = r_s+i;
     714        5562 :     ierr = DSVectors(d->eps->ds,DS_MAT_Y,&k,NULL);CHKERRQ(ierr);
     715        5562 :     ierr = DSGetArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
     716        5562 :     ierr = DSGetArray(d->eps->ds,DS_MAT_Y,&pY);CHKERRQ(ierr);
     717        5562 :     ierr = dvd_improvex_jd_proj_cuv(d,r_s+i,r_s+i+s,kr,data->theta,data->thetai,pX,pY,ld);CHKERRQ(ierr);
     718        5562 :     ierr = DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
     719        5562 :     ierr = DSRestoreArray(d->eps->ds,DS_MAT_Y,&pY);CHKERRQ(ierr);
     720             : 
     721             :     /* Check if the first eigenpairs are converged */
     722        5562 :     if (i == 0) {
     723        5025 :       PetscInt oldnpreconv = d->npreconv;
     724        5025 :       ierr = d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv);CHKERRQ(ierr);
     725        5025 :       if (d->npreconv > oldnpreconv) break;
     726             :     }
     727             : 
     728             :     /* Test the odd situation of solving Ax=b with A=I */
     729             : #if !defined(PETSC_USE_COMPLEX)
     730        5318 :     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        5318 :     if (data->ksp && !odd_situation) {
     736             :       /* kr <- -kr */
     737         470 :       for (j=0;j<s;j++) {
     738         470 :         ierr = VecScale(kr[j],-1.0);CHKERRQ(ierr);
     739             :       }
     740             : 
     741             :       /* Compose kr and D */
     742         470 :       kr0[0] = kr[0];
     743         470 :       kr0[1] = (s==2 ? kr[1] : NULL);
     744         470 :       ierr = VecCreateCompWithVecs(kr0,data->ksp_max_size,data->friends,&kr_comp);CHKERRQ(ierr);
     745         470 :       ierr = BVGetColumn(d->eps->V,kV+i,&D[0]);CHKERRQ(ierr);
     746         470 :       if (s==2) { ierr = BVGetColumn(d->eps->V,kV+i+1,&D[1]);CHKERRQ(ierr); }
     747         470 :       else D[1] = NULL;
     748         470 :       ierr = VecCreateCompWithVecs(D,data->ksp_max_size,data->friends,&D_comp);CHKERRQ(ierr);
     749         470 :       ierr = VecCompSetSubVecs(data->friends,s,NULL);CHKERRQ(ierr);
     750             : 
     751             :       /* Solve the correction equation */
     752         470 :       ierr = KSPSetTolerances(data->ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,maxits);CHKERRQ(ierr);
     753         470 :       ierr = KSPSolve(data->ksp,kr_comp,D_comp);CHKERRQ(ierr);
     754         470 :       ierr = KSPGetIterationNumber(data->ksp,&lits);CHKERRQ(ierr);
     755             : 
     756             :       /* Destroy the composed ks and D */
     757         470 :       ierr = VecDestroy(&kr_comp);CHKERRQ(ierr);
     758         470 :       ierr = VecDestroy(&D_comp);CHKERRQ(ierr);
     759         470 :       ierr = BVRestoreColumn(d->eps->V,kV+i,&D[0]);CHKERRQ(ierr);
     760         470 :       if (s==2) { ierr = BVRestoreColumn(d->eps->V,kV+i+1,&D[1]);CHKERRQ(ierr); }
     761             : 
     762             :     /* If GD */
     763             :     } else {
     764        4848 :       ierr = BVGetColumn(d->eps->V,kV+i,&D[0]);CHKERRQ(ierr);
     765        4848 :       if (s==2) { ierr = BVGetColumn(d->eps->V,kV+i+1,&D[1]);CHKERRQ(ierr); }
     766        4905 :       for (j=0;j<s;j++) {
     767        4905 :         ierr = d->improvex_precond(d,r_s+i+j,kr[j],D[j]);CHKERRQ(ierr);
     768             :       }
     769        4848 :       ierr = dvd_improvex_apply_proj(d,D,s);CHKERRQ(ierr);
     770        4848 :       ierr = BVRestoreColumn(d->eps->V,kV+i,&D[0]);CHKERRQ(ierr);
     771        4848 :       if (s==2) { ierr = BVRestoreColumn(d->eps->V,kV+i+1,&D[1]);CHKERRQ(ierr); }
     772             :     }
     773             :     /* Prevent that short vectors are discarded in the orthogonalization */
     774        5318 :     if (i == 0 && d->eps->errest[d->nconv+r_s] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s] < PETSC_MAX_REAL) {
     775        4832 :       for (j=0;j<s;j++) {
     776        4832 :         ierr = BVScaleColumn(d->eps->V,kV+i+j,1.0/d->eps->errest[d->nconv+r_s]);CHKERRQ(ierr);
     777             :       }
     778             :     }
     779        5318 :     ierr = SlepcVecPoolRestoreVecs(d->auxV,s,&kr);CHKERRQ(ierr);
     780             :   }
     781        5025 :   *size_D = i;
     782        5025 :   if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);
     783        5025 :   PetscFunctionReturn(0);
     784             : }
     785             : 
     786         219 : PetscErrorCode dvd_improvex_jd(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs,PetscBool dynamic)
     787             : {
     788             :   PetscErrorCode ierr;
     789             :   dvdImprovex_jd *data;
     790             :   PetscBool      useGD;
     791             :   PC             pc;
     792             :   PetscInt       size_P;
     793             : 
     794         219 :   PetscFunctionBegin;
     795             :   /* Setting configuration constrains */
     796         219 :   ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&useGD);CHKERRQ(ierr);
     797             : 
     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         219 :   if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) {
     802          60 :     max_bs++;
     803          60 :     b->max_size_P = PetscMax(b->max_size_P,2);
     804             :   } else
     805             : #endif
     806             :   {
     807         159 :     b->max_size_P = PetscMax(b->max_size_P,1);
     808             :   }
     809         219 :   b->max_size_X = PetscMax(b->max_size_X,max_bs);
     810         219 :   size_P = b->max_size_P;
     811             : 
     812             :   /* Setup the preconditioner */
     813         219 :   if (ksp) {
     814         219 :     ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
     815         219 :     ierr = dvd_static_precond_PC(d,b,pc);CHKERRQ(ierr);
     816             :   } else {
     817           0 :     ierr = dvd_static_precond_PC(d,b,0);CHKERRQ(ierr);
     818             :   }
     819             : 
     820             :   /* Setup the step */
     821         219 :   if (b->state >= DVD_STATE_CONF) {
     822          73 :     ierr = PetscNewLog(d->eps,&data);CHKERRQ(ierr);
     823          73 :     data->dynamic = dynamic;
     824          73 :     ierr = PetscMalloc1(size_P*size_P,&data->XKZ);CHKERRQ(ierr);
     825          73 :     ierr = PetscMalloc1(size_P*size_P,&data->iXKZ);CHKERRQ(ierr);
     826          73 :     ierr = PetscMalloc1(size_P,&data->iXKZPivots);CHKERRQ(ierr);
     827          73 :     data->ldXKZ = size_P;
     828          73 :     data->size_X = b->max_size_X;
     829          73 :     d->improveX_data = data;
     830          73 :     data->ksp = useGD? NULL: ksp;
     831          73 :     data->d = d;
     832          73 :     d->improveX = dvd_improvex_jd_gen;
     833             : #if !defined(PETSC_USE_COMPLEX)
     834          73 :     if (!DVD_IS(d->sEP,DVD_EP_HERMITIAN)) data->ksp_max_size = 2;
     835             :     else
     836             : #endif
     837          53 :       data->ksp_max_size = 1;
     838             :     /* Create various vector basis */
     839          73 :     ierr = BVDuplicateResize(d->eps->V,size_P,&data->KZ);CHKERRQ(ierr);
     840          73 :     ierr = BVSetMatrix(data->KZ,NULL,PETSC_FALSE);CHKERRQ(ierr);
     841          73 :     ierr = BVDuplicate(data->KZ,&data->U);CHKERRQ(ierr);
     842             : 
     843          73 :     ierr = EPSDavidsonFLAdd(&d->startList,dvd_improvex_jd_start);CHKERRQ(ierr);
     844          73 :     ierr = EPSDavidsonFLAdd(&d->endList,dvd_improvex_jd_end);CHKERRQ(ierr);
     845          73 :     ierr = EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_jd_d);CHKERRQ(ierr);
     846             :   }
     847         219 :   PetscFunctionReturn(0);
     848             : }
     849             : 
     850             : #if !defined(PETSC_USE_COMPLEX)
     851          59 : 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             : {
     853             :   PetscErrorCode ierr;
     854             :   PetscScalar    rAr,iAr,rAi,iAi,rBr,iBr,rBi,iBi,b0,b2,b4,b6,b7;
     855             : 
     856          59 :   PetscFunctionBegin;
     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          59 :   ierr = VecDotBegin(Axr,ur,&rAr);CHKERRQ(ierr); /* r*A*r */
     861          59 :   ierr = VecDotBegin(Axr,ui,&iAr);CHKERRQ(ierr); /* i*A*r */
     862          59 :   ierr = VecDotBegin(Axi,ur,&rAi);CHKERRQ(ierr); /* r*A*i */
     863          59 :   ierr = VecDotBegin(Axi,ui,&iAi);CHKERRQ(ierr); /* i*A*i */
     864          59 :   ierr = VecDotBegin(Bxr,ur,&rBr);CHKERRQ(ierr); /* r*B*r */
     865          59 :   ierr = VecDotBegin(Bxr,ui,&iBr);CHKERRQ(ierr); /* i*B*r */
     866          59 :   ierr = VecDotBegin(Bxi,ur,&rBi);CHKERRQ(ierr); /* r*B*i */
     867          59 :   ierr = VecDotBegin(Bxi,ui,&iBi);CHKERRQ(ierr); /* i*B*i */
     868          59 :   ierr = VecDotEnd(Axr,ur,&rAr);CHKERRQ(ierr); /* r*A*r */
     869          59 :   ierr = VecDotEnd(Axr,ui,&iAr);CHKERRQ(ierr); /* i*A*r */
     870          59 :   ierr = VecDotEnd(Axi,ur,&rAi);CHKERRQ(ierr); /* r*A*i */
     871          59 :   ierr = VecDotEnd(Axi,ui,&iAi);CHKERRQ(ierr); /* i*A*i */
     872          59 :   ierr = VecDotEnd(Bxr,ur,&rBr);CHKERRQ(ierr); /* r*B*r */
     873          59 :   ierr = VecDotEnd(Bxr,ui,&iBr);CHKERRQ(ierr); /* i*B*r */
     874          59 :   ierr = VecDotEnd(Bxi,ur,&rBi);CHKERRQ(ierr); /* r*B*i */
     875          59 :   ierr = VecDotEnd(Bxi,ui,&iBi);CHKERRQ(ierr); /* i*B*i */
     876          59 :   b0 = rAr+iAi; /* rAr+iAi */
     877          59 :   b2 = rAi-iAr; /* rAi-iAr */
     878          59 :   b4 = rBr+iBi; /* rBr+iBi */
     879          59 :   b6 = rBi-iBr; /* rBi-iBr */
     880          59 :   b7 = b4*b4 + b6*b6; /* k */
     881          59 :   *eigr = (b0*b4 + b2*b6) / b7; /* eig_r */
     882          59 :   *eigi = (b2*b4 - b0*b6) / b7; /* eig_i */
     883          59 :   PetscFunctionReturn(0);
     884             : }
     885             : #endif
     886             : 
     887        5562 : PETSC_STATIC_INLINE PetscErrorCode dvd_compute_n_rr(PetscInt i_s,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,Vec *u,Vec *Ax,Vec *Bx)
     888             : {
     889             :   PetscErrorCode ierr;
     890             :   PetscInt       i;
     891             :   PetscScalar    b0,b1;
     892             : 
     893        5562 :   PetscFunctionBegin;
     894        5562 :   for (i=0; i<n; i++) {
     895             : #if !defined(PETSC_USE_COMPLEX)
     896        5562 :     if (eigi[i_s+i] != 0.0) {
     897          59 :       PetscScalar eigr0=0.0,eigi0=0.0;
     898          59 :       ierr = dvd_complex_rayleigh_quotient(u[i],u[i+1],Ax[i],Ax[i+1],Bx[i],Bx[i+1],&eigr0,&eigi0);CHKERRQ(ierr);
     899          59 :       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           0 :         ierr = 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);CHKERRQ(ierr);
     901             :       }
     902          59 :       i++;
     903             :     } else
     904             : #endif
     905             :     {
     906        5503 :       ierr = VecDotBegin(Ax[i],u[i],&b0);CHKERRQ(ierr);
     907        5503 :       ierr = VecDotBegin(Bx[i],u[i],&b1);CHKERRQ(ierr);
     908        5503 :       ierr = VecDotEnd(Ax[i],u[i],&b0);CHKERRQ(ierr);
     909        5503 :       ierr = VecDotEnd(Bx[i],u[i],&b1);CHKERRQ(ierr);
     910        5503 :       b0 = b0/b1;
     911        5503 :       if (PetscAbsScalar(eigr[i_s+i]-b0)/PetscAbsScalar(eigr[i_s+i]) > 1e-10) {
     912           0 :         ierr = 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));CHKERRQ(ierr);
     913             :       }
     914             :     }
     915             :   }
     916        5562 :   PetscFunctionReturn(0);
     917             : }
     918             : 
     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        5562 : 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             : {
     928             :   PetscErrorCode ierr;
     929        5562 :   PetscInt       n = i_e-i_s,i;
     930             :   PetscScalar    *b;
     931             :   Vec            *Ax,*Bx,*r;
     932             :   Mat            M;
     933             :   BV             X;
     934             : 
     935        5562 :   PetscFunctionBegin;
     936        5562 :   ierr = BVDuplicateResize(d->eps->V,4,&X);CHKERRQ(ierr);
     937        5562 :   ierr = MatCreateSeqDense(PETSC_COMM_SELF,4,4,NULL,&M);CHKERRQ(ierr);
     938             :   /* u <- X(i) */
     939        5562 :   ierr = dvd_improvex_compute_X(d,i_s,i_e,u,pX,ld);CHKERRQ(ierr);
     940             : 
     941             :   /* v <- theta[0]A*u + theta[1]*B*u */
     942             : 
     943             :   /* Bx <- B*X(i) */
     944        5562 :   Bx = kr;
     945        5562 :   if (d->BX) {
     946        1419 :     for (i=i_s; i<i_e; ++i) {
     947        1419 :       ierr = BVMultVec(d->BX,1.0,0.0,Bx[i-i_s],&pX[ld*i]);CHKERRQ(ierr);
     948             :     }
     949             :   } else {
     950        4202 :     for (i=0;i<n;i++) {
     951        4202 :       if (d->B) {
     952           0 :         ierr = MatMult(d->B, u[i], Bx[i]);CHKERRQ(ierr);
     953             :       } else {
     954        4202 :         ierr = VecCopy(u[i], Bx[i]);CHKERRQ(ierr);
     955             :       }
     956             :     }
     957             :   }
     958             : 
     959             :   /* Ax <- A*X(i) */
     960        5562 :   ierr = SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&r);CHKERRQ(ierr);
     961        5562 :   Ax = r;
     962       11183 :   for (i=i_s; i<i_e; ++i) {
     963        5621 :     ierr = BVMultVec(d->AX,1.0,0.0,Ax[i-i_s],&pX[ld*i]);CHKERRQ(ierr);
     964             :   }
     965             : 
     966             :   /* v <- Y(i) */
     967        5621 :   for (i=i_s; i<i_e; ++i) {
     968        5621 :     ierr = BVMultVec(d->W?d->W:d->eps->V,1.0,0.0,v[i-i_s],&pY[ld*i]);CHKERRQ(ierr);
     969             :   }
     970             : 
     971             :   /* Recompute the eigenvalue */
     972        5562 :   ierr = dvd_compute_n_rr(i_s,n,d->eigr,d->eigi,v,Ax,Bx);CHKERRQ(ierr);
     973             : 
     974        5562 :   for (i=0;i<n;i++) {
     975             : #if !defined(PETSC_USE_COMPLEX)
     976        5562 :     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          59 :       ierr = MatDenseGetArray(M,&b);CHKERRQ(ierr);
     982          59 :       b[0] = b[5] = PetscConj(theta[2*i]);
     983          59 :       b[2] = b[7] = -theta[2*i+1];
     984          59 :       b[6] = -(b[3] = thetai[i]);
     985          59 :       b[1] = b[4] = 0.0;
     986          59 :       b[8] = b[13] = 1.0/d->nX[i_s+i];
     987          59 :       b[10] = b[15] = -d->eigr[i_s+i]/d->nX[i_s+i];
     988          59 :       b[14] = -(b[11] = d->eigi[i_s+i]/d->nX[i_s+i]);
     989          59 :       b[9] = b[12] = 0.0;
     990          59 :       ierr = MatDenseRestoreArray(M,&b);CHKERRQ(ierr);
     991          59 :       ierr = BVInsertVec(X,0,Ax[i]);CHKERRQ(ierr);
     992          59 :       ierr = BVInsertVec(X,1,Ax[i+1]);CHKERRQ(ierr);
     993          59 :       ierr = BVInsertVec(X,2,Bx[i]);CHKERRQ(ierr);
     994          59 :       ierr = BVInsertVec(X,3,Bx[i+1]);CHKERRQ(ierr);
     995          59 :       ierr = BVSetActiveColumns(X,0,4);CHKERRQ(ierr);
     996          59 :       ierr = BVMultInPlace(X,M,0,4);CHKERRQ(ierr);
     997          59 :       ierr = BVCopyVec(X,0,Ax[i]);CHKERRQ(ierr);
     998          59 :       ierr = BVCopyVec(X,1,Ax[i+1]);CHKERRQ(ierr);
     999          59 :       ierr = BVCopyVec(X,2,Bx[i]);CHKERRQ(ierr);
    1000          59 :       ierr = BVCopyVec(X,3,Bx[i+1]);CHKERRQ(ierr);
    1001          59 :       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        5503 :       ierr = MatDenseGetArray(M,&b);CHKERRQ(ierr);
    1008        5503 :       b[0] = PetscConj(theta[i*2]);
    1009        5503 :       b[1] = theta[i*2+1];
    1010        5503 :       b[4] = 1.0/d->nX[i_s+i];
    1011        5503 :       b[5] = -d->eigr[i_s+i]/d->nX[i_s+i];
    1012        5503 :       ierr = MatDenseRestoreArray(M,&b);CHKERRQ(ierr);
    1013        5503 :       ierr = BVInsertVec(X,0,Ax[i]);CHKERRQ(ierr);
    1014        5503 :       ierr = BVInsertVec(X,1,Bx[i]);CHKERRQ(ierr);
    1015        5503 :       ierr = BVSetActiveColumns(X,0,2);CHKERRQ(ierr);
    1016        5503 :       ierr = BVMultInPlace(X,M,0,2);CHKERRQ(ierr);
    1017        5503 :       ierr = BVCopyVec(X,0,Ax[i]);CHKERRQ(ierr);
    1018        5503 :       ierr = BVCopyVec(X,1,Bx[i]);CHKERRQ(ierr);
    1019             :     }
    1020             :   }
    1021        5621 :   for (i=0; i<n; i++) d->nX[i_s+i] = 1.0;
    1022             : 
    1023             :   /* v <- K^{-1} r = K^{-1}(theta_2i'*Ax + theta_2i+1*Bx) */
    1024        5621 :   for (i=0;i<n;i++) {
    1025        5621 :     ierr = d->improvex_precond(d,i_s+i,r[i],v[i]);CHKERRQ(ierr);
    1026             :   }
    1027        5562 :   ierr = SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&r);CHKERRQ(ierr);
    1028             : 
    1029             :   /* kr <- P*(Ax - eig_i*Bx) */
    1030        5562 :   ierr = d->calcpairs_proj_res(d,i_s,i_e,kr);CHKERRQ(ierr);
    1031        5562 :   ierr = BVDestroy(&X);CHKERRQ(ierr);
    1032        5562 :   ierr = MatDestroy(&M);CHKERRQ(ierr);
    1033        5562 :   PetscFunctionReturn(0);
    1034             : }
    1035             : 
    1036        5621 : static PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d,PetscInt i,PetscScalar* theta,PetscScalar* thetai,PetscInt *maxits,PetscReal *tol)
    1037             : {
    1038        5621 :   dvdImprovex_jd *data = (dvdImprovex_jd*)d->improveX_data;
    1039             :   PetscReal      a;
    1040             : 
    1041        5621 :   PetscFunctionBegin;
    1042        5621 :   a = SlepcAbsEigenvalue(d->eigr[i],d->eigi[i]);
    1043             : 
    1044        5621 :   if (d->nR[i] < data->fix*a) {
    1045         382 :     theta[0] = d->eigr[i];
    1046         382 :     theta[1] = 1.0;
    1047             : #if !defined(PETSC_USE_COMPLEX)
    1048         382 :     *thetai = d->eigi[i];
    1049             : #endif
    1050             :   } else {
    1051        5239 :     theta[0] = d->target[0];
    1052        5239 :     theta[1] = d->target[1];
    1053             : #if !defined(PETSC_USE_COMPLEX)
    1054        5239 :     *thetai = 0.0;
    1055             : #endif
    1056             : }
    1057             : 
    1058             : #if defined(PETSC_USE_COMPLEX)
    1059             :   if (thetai) *thetai = 0.0;
    1060             : #endif
    1061        5621 :   *maxits = data->maxits;
    1062        5621 :   *tol = data->tol;
    1063        5621 :   PetscFunctionReturn(0);
    1064             : }
    1065             : 
    1066         219 : PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard *d,dvdBlackboard *b,PetscInt maxits,PetscReal tol,PetscReal fix)
    1067             : {
    1068         219 :   dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;
    1069             : 
    1070         219 :   PetscFunctionBegin;
    1071             :   /* Setup the step */
    1072         219 :   if (b->state >= DVD_STATE_CONF) {
    1073          73 :     data->maxits = maxits;
    1074          73 :     data->tol = tol;
    1075          73 :     data->fix = fix;
    1076          73 :     d->improvex_jd_lit = dvd_improvex_jd_lit_const_0;
    1077             :   }
    1078         219 :   PetscFunctionReturn(0);
    1079             : }
    1080             : 
    1081         219 : PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d,dvdBlackboard *b)
    1082             : {
    1083         219 :   PetscFunctionBegin;
    1084             :   /* Setup the step */
    1085         219 :   if (b->state >= DVD_STATE_CONF) {
    1086          73 :     d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX;
    1087             :   }
    1088         219 :   PetscFunctionReturn(0);
    1089             : }
    1090             : 
    1091       24499 : PetscErrorCode dvd_improvex_compute_X(dvdDashboard *d,PetscInt i_s,PetscInt i_e,Vec *u_,PetscScalar *pX,PetscInt ld)
    1092             : {
    1093             :   PetscErrorCode ierr;
    1094       24499 :   PetscInt       n = i_e - i_s,i;
    1095             :   Vec            *u;
    1096             : 
    1097       24499 :   PetscFunctionBegin;
    1098       24499 :   if (u_) u = u_;
    1099        2778 :   else if (d->correctXnorm) {
    1100         183 :     ierr = SlepcVecPoolGetVecs(d->auxV,i_e-i_s,&u);CHKERRQ(ierr);
    1101             :   }
    1102       24499 :   if (u_ || d->correctXnorm) {
    1103       21963 :     for (i=0; i<n; i++) {
    1104       21963 :       ierr = BVMultVec(d->eps->V,1.0,0.0,u[i],&pX[ld*(i+i_s)]);CHKERRQ(ierr);
    1105             :     }
    1106             :   }
    1107             :   /* nX(i) <- ||X(i)|| */
    1108       24499 :   if (d->correctXnorm) {
    1109        1520 :     for (i=0; i<n; i++) {
    1110        1520 :       ierr = VecNormBegin(u[i],NORM_2,&d->nX[i_s+i]);CHKERRQ(ierr);
    1111             :     }
    1112        1520 :     for (i=0; i<n; i++) {
    1113        1520 :       ierr = VecNormEnd(u[i],NORM_2,&d->nX[i_s+i]);CHKERRQ(ierr);
    1114             :     }
    1115             : #if !defined(PETSC_USE_COMPLEX)
    1116        1520 :     for (i=0;i<n;i++) {
    1117        1520 :       if (d->eigi[i_s+i] != 0.0) {
    1118           0 :         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           0 :         i++;
    1120             :       }
    1121             :     }
    1122             : #endif
    1123             :   } else {
    1124       23065 :     for (i=0;i<n;i++) d->nX[i_s+i] = 1.0;
    1125             :   }
    1126       24499 :   if (d->correctXnorm && !u_) {
    1127         183 :     ierr = SlepcVecPoolRestoreVecs(d->auxV,i_e-i_s,&u);CHKERRQ(ierr);
    1128             :   }
    1129       24499 :   PetscFunctionReturn(0);
    1130             : }
    1131             : 

Generated by: LCOV version 1.13