Subversion Repositories slepc-dev

/*

  SLEPc eigensolver: "davidson"

  Step: improve the eigenvectors X

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

#include "davidson.h"

#include <slepcvec.h>

#include <slepcblaslapack.h>

PetscErrorCode dvd_improvex_PfuncV(dvdDashboard *d, void *funcV, Vec *D,

  PetscInt max_size_D, PetscInt r_s, PetscInt r_e, Vec *auxV,

  PetscScalar *auxS);

PetscErrorCode dvd_matmult_jd(Mat A, Vec in, Vec out);

PetscErrorCode dvd_matgetvecs_jd(Mat A, Vec *right, Vec *left);

PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d);

PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d);

PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d);

PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d, Vec *D,

                                   PetscInt max_size_D, PetscInt r_s,

                                   PetscInt r_e, PetscInt *size_D);

PetscErrorCode dvd_improvex_jd_proj_cuv(dvdDashboard *d, PetscInt i_s,

  PetscInt i_e, Vec **u, Vec **v, Vec *kr, Vec **auxV, PetscScalar **auxS,

  PetscScalar *theta, PetscScalar *thetai, PetscScalar *pX, PetscScalar *pY,

  PetscInt ld);

PetscErrorCode dvd_improvex_jd_proj_uv_KXX(dvdDashboard *d, PetscInt i_s,

  PetscInt i_e, Vec *u, Vec *v, Vec *kr, Vec *auxV, PetscScalar *theta,

  PetscScalar *thetai, PetscScalar *pX, PetscScalar *pY, PetscInt ld);

PetscErrorCode dvd_improvex_jd_proj_uv_KZX(dvdDashboard *d, PetscInt i_s,

  PetscInt i_e, Vec *u, Vec *v, Vec *kr, Vec *auxV, PetscScalar *theta,

  PetscScalar *thetai, PetscScalar *pX, PetscScalar *pY, PetscInt ld);

PetscErrorCode dvd_improvex_jd_lit_const_0(dvdDashboard *d, PetscInt i,

  PetscScalar* theta, PetscScalar* thetai,

  PetscInt *maxits, PetscReal *tol);

PetscErrorCode dvd_improvex_get_eigenvectors(dvdDashboard *d, PetscScalar *pX,

  PetscScalar *pY, PetscInt ld_,

  PetscScalar *auxS, PetscInt size_auxS);

PetscErrorCode dvd_improvex_apply_proj(dvdDashboard *d, Vec *V, PetscInt cV, PetscScalar *auxS);

#define size_Z (64*4)

/**** JD update step (I - Kfg'/(g'Kf)) K(A - sB) (I - Kfg'/(g'Kf)) t = (I - Kfg'/(g'Kf))r  ****/

typedef struct {

  PetscInt size_X;

  void

    *old_improveX_data;   /* old improveX_data */

  improveX_type

    old_improveX;         /* old improveX */

  KSP ksp;                /* correction equation solver */

  Vec

    friends,              /* reference vector for composite vectors */

    *auxV;                /* auxiliar vectors */

  PetscScalar *auxS,      /* auxiliar scalars */

    *theta,

    *thetai;              /* the shifts used in the correction eq. */

  PetscInt maxits,        /* maximum number of iterations */

    r_s, r_e,             /* the selected eigenpairs to improve */

    ksp_max_size;         /* the ksp maximum subvectors size */

  PetscReal tol,          /* the maximum solution tolerance */

    lastTol,              /* last tol for dynamic stopping criterion */

    fix;                  /* tolerance for using the approx. eigenvalue */

  PetscBool

    dynamic;              /* if the dynamic stopping criterion is applied */

  dvdDashboard

    *d;                   /* the currect dvdDashboard reference */

  PC old_pc;              /* old pc in ksp */

  Vec

    *u,                   /* new X vectors */

    *real_KZ,             /* original KZ */

    *KZ;                  /* KZ vecs for the projector KZ*inv(X'*KZ)*X' */

  PetscScalar

   *XKZ,                  /* X'*KZ */

   *iXKZ;                 /* inverse of XKZ */

  PetscInt

    ldXKZ,                /* leading dimension of XKZ */

    size_iXKZ,            /* size of iXKZ */

    ldiXKZ,               /* leading dimension of iXKZ */

    size_KZ,              /* size of converged KZ */

    size_real_KZ,         /* original size of KZ */

    size_cX,              /* last value of d->size_cX */

    old_size_X;           /* last number of improved vectors */

  PetscBLASInt

    *iXKZPivots;          /* array of pivots */

} dvdImprovex_jd;

#define _Ceil(A,B) ((A)/(B)+((A)%(B)==0?0:1))

#define FromIntToScalar(S) (_Ceil((S)*sizeof(PetscBLASInt),sizeof(PetscScalar)))

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd"

PetscErrorCode dvd_improvex_jd(dvdDashboard *d, dvdBlackboard *b, KSP ksp,

                               PetscInt max_bs, PetscInt cX_impr, PetscBool dynamic)

{

  PetscErrorCode  ierr;

  dvdImprovex_jd  *data;

  PetscBool       t, herm = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;

  PC              pc;

  PetscInt        size_P;

  PetscFunctionBegin;

  /* Setting configuration constrains */

  ierr = PetscTypeCompare((PetscObject)ksp, KSPPREONLY, &t); CHKERRQ(ierr);

  if (t) ksp = PETSC_NULL;

  /* If the arithmetic is real and the problem is not Hermitian, then

     the block size is incremented in one */

#if !defined(PETSC_USE_COMPLEX)

  if (!herm) {

    max_bs++;

    b->max_size_P = PetscMax(b->max_size_P, 2);

  } else

#endif

    b->max_size_P = PetscMax(b->max_size_P, 1);

  b->max_size_X = PetscMax(b->max_size_X, max_bs);

  size_P = b->max_size_P+cX_impr;

  b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X*(ksp?3:2));

                                                            /* u, kr?, auxV */

  b->own_scalars+= size_P*size_P; /* XKZ */

  b->max_size_auxS = PetscMax(b->max_size_auxS,

    (herm?0:1)*2*b->max_size_proj*b->max_size_proj + /* pX, pY */

    b->max_size_X*3 + /* theta, thetai */

    size_P*size_P + /* iXKZ */

    FromIntToScalar(size_P) + /* iXkZPivots */

    PetscMax(PetscMax(

      3*b->max_size_proj*b->max_size_X, /* dvd_improvex_apply_proj */

      8*cX_impr*b->max_size_X),

                                          /* dvd_improvex_jd_proj_cuv_KZX */

      (herm?0:1)*11*b->max_size_proj+4*b->max_size_proj*b->max_size_proj));

                                           /* dvd_improvex_get_eigenvectors */

  b->own_vecs+= size_P; /* KZ */

  /* Setup the preconditioner */

  if (ksp) {

    ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);

    ierr = dvd_static_precond_PC(d, b, pc); CHKERRQ(ierr);

  } else {

    ierr = dvd_static_precond_PC(d, b, 0); CHKERRQ(ierr);

  }

  /* Setup the step */

  if (b->state >= DVD_STATE_CONF) {

    ierr = PetscMalloc(sizeof(dvdImprovex_jd), &data); CHKERRQ(ierr);

    data->dynamic = dynamic;

    data->size_real_KZ = size_P;

    data->real_KZ = b->free_vecs; b->free_vecs+= data->size_real_KZ;

    d->max_size_cX_in_impr = cX_impr;

    data->XKZ = b->free_scalars; b->free_scalars+= size_P*size_P;

    data->ldXKZ = size_P;

    data->size_X = b->max_size_X;

    data->old_improveX_data = d->improveX_data;

    d->improveX_data = data;

    data->old_improveX = d->improveX;

    data->ksp = ksp;

    data->d = d;

    d->improveX = dvd_improvex_jd_gen;

    data->ksp_max_size = max_bs;

    DVD_FL_ADD(d->startList, dvd_improvex_jd_start);

    DVD_FL_ADD(d->endList, dvd_improvex_jd_end);

    DVD_FL_ADD(d->destroyList, dvd_improvex_jd_d);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_start"

PetscErrorCode dvd_improvex_jd_start(dvdDashboard *d)

{

  PetscErrorCode  ierr;

  dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

  PetscInt        rA, cA, rlA, clA;

  Mat             A;

  PetscBool       t;

  PC              pc;

  PetscFunctionBegin;

  data->KZ = data->real_KZ;

  data->size_KZ = data->size_cX = data->old_size_X = 0;

  data->lastTol = data->dynamic?0.5:0.0;

  /* Setup the ksp */

  if(data->ksp) {

    /* Create the reference vector */

    ierr = VecCreateCompWithVecs(d->V, data->ksp_max_size, PETSC_NULL,

                                 &data->friends); CHKERRQ(ierr);

    /* Save the current pc and set a PCNONE */

    ierr = KSPGetPC(data->ksp, &data->old_pc); CHKERRQ(ierr);

    ierr = PetscTypeCompare((PetscObject)data->old_pc, PCNONE, &t);

    CHKERRQ(ierr);

    data->lastTol = 0.5;

    if (t) {

      data->old_pc = 0;

    } else {

      ierr = PetscObjectReference((PetscObject)data->old_pc); CHKERRQ(ierr);

      ierr = PCCreate(((PetscObject)d->eps)->comm, &pc); CHKERRQ(ierr);

      ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr);

      ierr = PCSetOperators(pc, d->A, d->A, SAME_PRECONDITIONER); CHKERRQ(ierr);

      ierr = KSPSetPC(data->ksp, pc); CHKERRQ(ierr);

      ierr = PCDestroy(&pc); CHKERRQ(ierr);

    }

    /* Create the (I-v*u')*K*(A-s*B) matrix */

    ierr = MatGetSize(d->A, &rA, &cA); CHKERRQ(ierr);

    ierr = MatGetLocalSize(d->A, &rlA, &clA); CHKERRQ(ierr);

    ierr = MatCreateShell(((PetscObject)d->A)->comm, rlA*data->ksp_max_size,

                          clA*data->ksp_max_size, rA*data->ksp_max_size,

                          cA*data->ksp_max_size, data, &A); CHKERRQ(ierr);

    ierr = MatShellSetOperation(A, MATOP_MULT,

                                (void(*)(void))dvd_matmult_jd); CHKERRQ(ierr);

    ierr = MatShellSetOperation(A, MATOP_GET_VECS,

                                (void(*)(void))dvd_matgetvecs_jd);

    CHKERRQ(ierr);

    /* Try to avoid KSPReset */

    ierr = KSPGetOperatorsSet(data->ksp,&t,PETSC_NULL);CHKERRQ(ierr);

    if (t) {

      Mat M;

      PetscInt rM;

      ierr = KSPGetOperators(data->ksp,&M,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);

      ierr = MatGetSize(M,&rM,PETSC_NULL);CHKERRQ(ierr);

      if (rM != rA*data->ksp_max_size) { ierr = KSPReset(data->ksp);CHKERRQ(ierr); }

    }

    ierr = KSPSetOperators(data->ksp, A, A, SAME_PRECONDITIONER);

    CHKERRQ(ierr);

    ierr = KSPSetUp(data->ksp); CHKERRQ(ierr);

    ierr = MatDestroy(&A); CHKERRQ(ierr);

  } else {

    data->old_pc = 0;

    data->friends = PETSC_NULL;

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_end"

PetscErrorCode dvd_improvex_jd_end(dvdDashboard *d)

{

  PetscErrorCode  ierr;

  dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

  PetscFunctionBegin;

  if (data->friends) { ierr = VecDestroy(&data->friends); CHKERRQ(ierr); }

  /* Restore the pc of ksp */

  if (data->old_pc) {

    ierr = KSPSetPC(data->ksp, data->old_pc); CHKERRQ(ierr);

    ierr = PCDestroy(&data->old_pc); CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_d"

PetscErrorCode dvd_improvex_jd_d(dvdDashboard *d)

{

  PetscErrorCode  ierr;

  dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

  PetscFunctionBegin;

  /* Restore changes in dvdDashboard */

  d->improveX_data = data->old_improveX_data;

  /* Free local data and objects */

  ierr = PetscFree(data); CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_gen"

PetscErrorCode dvd_improvex_jd_gen(dvdDashboard *d, Vec *D,

                             PetscInt max_size_D, PetscInt r_s,

                             PetscInt r_e, PetscInt *size_D)

{

  dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

  PetscErrorCode  ierr;

  PetscInt        i, j, n, maxits, maxits0, lits, s;

  PetscScalar     *pX, *pY, *auxS = d->auxS, *auxS0;

  PetscReal       tol, tol0;

  Vec             *u, *v, *kr, kr_comp, D_comp;

  PetscFunctionBegin;

  /* Quick exit */

  if ((max_size_D == 0) || r_e-r_s <= 0) {

   /* Callback old improveX */

    if (data->old_improveX) {

      d->improveX_data = data->old_improveX_data;

      data->old_improveX(d, PETSC_NULL, 0, 0, 0, PETSC_NULL);

      d->improveX_data = data;

    }

    PetscFunctionReturn(0);

  }

  n = PetscMin(PetscMin(data->size_X, max_size_D), r_e-r_s);

  if (n == 0) SETERRQ(PETSC_COMM_SELF,1, "n == 0!\n");

  if (data->size_X < r_e-r_s) SETERRQ(PETSC_COMM_SELF,1, "size_X < r_e-r_s!\n");

  /* Compute the eigenvectors of the selected pairs */

  if (DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {

    pX = pY = d->pX;

  } else {

    pX = auxS; auxS+= d->size_H*d->size_H;

    pY = auxS; auxS+= d->size_H*d->size_H;

    ierr = dvd_improvex_get_eigenvectors(d, pX, pY, d->size_H, auxS,

                                         d->size_auxS-(auxS-d->auxS));

    CHKERRQ(ierr);

  }

  /* Restart lastTol if a new pair converged */

  if (data->dynamic && data->size_cX < d->size_cX)

    data->lastTol = 0.5;

  for(i=0, s=0, auxS0=auxS; i<n; i+=s) {

    /* If the selected eigenvalue is complex, but the arithmetic is real... */

#if !defined(PETSC_USE_COMPLEX)

    if (PetscAbsScalar(d->eigi[i] != 0.0)) {

      if (i+2 <= max_size_D) s=2; else break;

    } else

#endif

      s=1;

    data->auxV = d->auxV;

    data->r_s = r_s+i; data->r_e = r_s+i+s;

    auxS = auxS0;

    data->theta = auxS; auxS+= 2*s;

    data->thetai = auxS; auxS+= s;

    /* If GD, kr = D */

    if (!data->ksp) {

      kr = &D[i];

    /* If JD, kr = auxV */

    } else {

      kr = data->auxV; data->auxV+= s;

    }

    /* Compute theta, maximum iterations and tolerance */

    maxits = 0; tol = 1;

    for(j=0; j<s; j++) {

      ierr = d->improvex_jd_lit(d, r_s+i+j, &data->theta[2*j],

                                &data->thetai[j], &maxits0, &tol0);

      CHKERRQ(ierr);

      maxits+= maxits0; tol*= tol0;

    }

    maxits/= s; tol = data->dynamic?data->lastTol:exp(log(tol)/s);

    /* Compute u, v and kr */

    ierr = dvd_improvex_jd_proj_cuv(d, r_s+i, r_s+i+s, &u, &v, kr,

             &data->auxV, &auxS, data->theta, data->thetai,

             &pX[d->size_H*(r_s+i+d->cX_in_H)], &pY[d->size_H*(r_s+i+d->cX_in_H)], d->size_H);

    CHKERRQ(ierr);

    data->u = u;

    /* Check if the first eigenpairs are converged */

    if (i == 0) {

      ierr = d->preTestConv(d, 0, s, s, PETSC_NULL, PETSC_NULL, &d->npreconv);

      CHKERRQ(ierr);

      if (d->npreconv > 0) break;

    }

    /* Compute kr <- kr - v*(u'*kr) */

    ierr = dvd_improvex_apply_proj(d, kr, s, auxS); CHKERRQ(ierr);

    /* If JD */

    if (data->ksp) {

      data->auxS = auxS;

      /* kr <- -kr */

      for(j=0; j<s; j++) {

        ierr = VecScale(kr[j], -1.0); CHKERRQ(ierr);

      }

      /* Compouse kr and D */

      ierr = VecCreateCompWithVecs(kr, data->ksp_max_size, data->friends,

                                   &kr_comp); CHKERRQ(ierr);

      ierr = VecCreateCompWithVecs(&D[i], data->ksp_max_size, data->friends,

                                   &D_comp); CHKERRQ(ierr);

      ierr = VecCompSetVecs(data->friends, PETSC_NULL, s); CHKERRQ(ierr);

      /* Solve the correction equation */

      ierr = KSPSetTolerances(data->ksp, tol, PETSC_DEFAULT, PETSC_DEFAULT,

                              maxits); CHKERRQ(ierr);

      ierr = KSPSolve(data->ksp, kr_comp, D_comp); CHKERRQ(ierr);

      ierr = KSPGetIterationNumber(data->ksp, &lits); CHKERRQ(ierr);

      d->eps->OP->lineariterations+= lits;

      /* Destroy the composed ks and D */

      ierr = VecDestroy(&kr_comp); CHKERRQ(ierr);

      ierr = VecDestroy(&D_comp); CHKERRQ(ierr);

    }

  }

  *size_D = i;

  if (data->dynamic) data->lastTol = PetscMax(data->lastTol/2.0,PETSC_MACHINE_EPSILON*10.0);

  /* Callback old improveX */

  if (data->old_improveX) {

    d->improveX_data = data->old_improveX_data;

    data->old_improveX(d, PETSC_NULL, 0, 0, 0, PETSC_NULL);

    d->improveX_data = data;

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_matmult_jd"

PetscErrorCode dvd_matmult_jd(Mat A, Vec in, Vec out)

{

  PetscErrorCode  ierr;

  dvdImprovex_jd  *data;

  PetscInt        n, i;

  const Vec       *inx, *outx, *Bx;

  PetscFunctionBegin;

  ierr = MatShellGetContext(A, (void**)&data); CHKERRQ(ierr);

  ierr = VecCompGetVecs(in, &inx, PETSC_NULL); CHKERRQ(ierr);

  ierr = VecCompGetVecs(out, &outx, PETSC_NULL); CHKERRQ(ierr);

  n = data->r_e - data->r_s;

  /* aux <- theta[1]A*in - theta[0]*B*in */

  for(i=0; i<n; i++) {

    ierr = MatMult(data->d->A, inx[i], data->auxV[i]); CHKERRQ(ierr);

  }

  if (data->d->B) {

    for(i=0; i<n; i++) {

      ierr = MatMult(data->d->B, inx[i], outx[i]); CHKERRQ(ierr);

    }

    Bx = outx;

  } else

    Bx = inx;

  for(i=0; i<n; i++) {

#if !defined(PETSC_USE_COMPLEX)

    if(data->d->eigi[data->r_s+i] != 0.0) {

      /* aux_i   <- [ t_2i+1*A*inx_i   - t_2i*Bx_i + ti_i*Bx_i+1;

         aux_i+1      t_2i+1*A*inx_i+1 - ti_i*Bx_i - t_2i*Bx_i+1  ] */

      ierr = VecAXPBYPCZ(data->auxV[i], -data->theta[2*i], data->thetai[i],

                         data->theta[2*i+1], Bx[i], Bx[i+1]);

      CHKERRQ(ierr);

      ierr = VecAXPBYPCZ(data->auxV[i+1], -data->thetai[i],

                         -data->theta[2*i], data->theta[2*i+1], Bx[i],

                         Bx[i+1]); CHKERRQ(ierr);

      i++;

    } else

#endif

    {

      ierr = VecAXPBY(data->auxV[i], -data->theta[i*2], data->theta[i*2+1],

                      Bx[i]); CHKERRQ(ierr);

    }

  }

  /* out <- K * aux */

  for(i=0; i<n; i++) {

    ierr = data->d->improvex_precond(data->d, data->r_s+i, data->auxV[i],

                                     outx[i]); CHKERRQ(ierr);

  }

  /* out <- out - v*(u'*out) */

  ierr = dvd_improvex_apply_proj(data->d, (Vec*)outx, n, data->auxS);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_matgetvecs_jd"

PetscErrorCode dvd_matgetvecs_jd(Mat A, Vec *right, Vec *left)

{

  PetscErrorCode  ierr;

  Vec             *r, *l;

  dvdImprovex_jd  *data;

  PetscInt        n, i;

  PetscFunctionBegin;

  ierr = MatShellGetContext(A, (void**)&data); CHKERRQ(ierr);

  n = data->ksp_max_size;

  if (right) {

    ierr = PetscMalloc(sizeof(Vec)*n, &r); CHKERRQ(ierr);

  }

  if (left) {

    ierr = PetscMalloc(sizeof(Vec)*n, &l); CHKERRQ(ierr);

  }

  for (i=0; i<n; i++) {

    ierr = MatGetVecs(data->d->A, right?&r[i]:PETSC_NULL,

                      left?&l[i]:PETSC_NULL); CHKERRQ(ierr);

  }

  if(right) {

    ierr = VecCreateCompWithVecs(r, n, data->friends, right); CHKERRQ(ierr);

    for (i=0; i<n; i++) {

      ierr = VecDestroy(&r[i]); CHKERRQ(ierr);

    }

  }

  if(left) {

    ierr = VecCreateCompWithVecs(l, n, data->friends, left); CHKERRQ(ierr);

    for (i=0; i<n; i++) {

      ierr = VecDestroy(&l[i]); CHKERRQ(ierr);

    }

  }

  if (right) {

    ierr = PetscFree(r); CHKERRQ(ierr);

  }

  if (left) {

    ierr = PetscFree(l); CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_proj_uv"

PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard *d, dvdBlackboard *b,

                                       ProjType_t p)

{

  PetscFunctionBegin;

  /* Setup the step */

  if (b->state >= DVD_STATE_CONF) {

    switch(p) {

    case DVD_PROJ_KXX:

      d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KXX; break;

    case DVD_PROJ_KZX:

      d->improvex_jd_proj_uv = dvd_improvex_jd_proj_uv_KZX; break;

    }

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "dvd_improvex_jd_proj_cuv"

/*

  Compute: u <- X, v <- K*(theta[0]*A+theta[1]*B)*X,

  kr <- K^{-1}*(A-eig*B)*X, being X <- V*pX[i_s..i_e-1], Y <- W*pY[i_s..i_e-1]

  where

  auxV, 4*(i_e-i_s) auxiliar global vectors

  pX,pY, the right and left eigenvectors of the projected system

  ld, the leading dimension of pX and pY

  auxS: max(8*bs*max_size_cX_in_proj,size_V*size_V)

*/

PetscErrorCode dvd_improvex_jd_proj_cuv(dvdDashboard *d, PetscInt i_s,

  PetscInt i_e, Vec **u, Vec **v, Vec *kr, Vec **auxV, PetscScalar **auxS,

  PetscScalar *theta, PetscScalar *thetai, PetscScalar *pX, PetscScalar *pY,

  PetscInt ld)

{

  PetscErrorCode  ierr;

  PetscInt        n = i_e - i_s, size_KZ, V_new, rm, i;

  PetscScalar     *wS0, *wS1;

  dvdImprovex_jd  *data = (dvdImprovex_jd*)d->improveX_data;

  PetscBLASInt    s, ldXKZ, info;

  PetscFunctionBegin;

  /* Check consistency */

  V_new = d->size_cX - data->size_cX;

  if (V_new > data->old_size_X) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }

  data->old_size_X = n;

  /* KZ <- KZ(rm:) */

  rm = PetscMax(V_new+data->size_KZ-d->max_size_cX_in_impr, 0);

  if (rm > 0) for (i=rm; i<data->size_KZ+V_new; i++) {

    ierr = VecCopy(data->KZ[i], data->KZ[i-rm]);CHKERRQ(ierr);

  }

  data->size_KZ+= V_new-rm;

  data->size_cX = d->size_cX;

  /* XKZ <- XKZ(rm:,rm:) */

  if (rm > 0) for (i=rm; i<data->size_KZ+V_new; i++) {

    ierr = SlepcDenseCopy(&data->XKZ[i*data->ldXKZ+i], data->ldXKZ, &data->XKZ[(i-rm)*data->ldXKZ+i-rm], data->ldXKZ, data->size_KZ, 1);CHKERRQ(ierr);

  }

  /* Compute X, KZ and KR */

  *u = *auxV; *auxV+= n;

  *v = &data->KZ[data->size_KZ];

  ierr = d->improvex_jd_proj_uv(d, i_s, i_e, *u, *v, kr, *auxV, theta, thetai,

                                pX, pY, ld);CHKERRQ(ierr);

  /* XKZ <- X'*KZ */

  size_KZ = data->size_KZ+n;

  wS0 = *auxS; wS1 = wS0+2*n*data->size_KZ+n*n;

  ierr = VecsMult(data->XKZ, 0, data->ldXKZ, d->V-data->size_KZ, 0, data->size_KZ, data->KZ, 0, size_KZ, wS0, wS1);CHKERRQ(ierr);

  ierr = VecsMult(&data->XKZ[data->size_KZ], 0, data->ldXKZ, *u, 0, n, data->KZ, 0, size_KZ, wS0, wS1);CHKERRQ(ierr);

  /* iXKZ <- inv(XKZ) */

  s = PetscBLASIntCast(size_KZ);

  data->ldiXKZ = data->size_iXKZ = size_KZ;

  data->iXKZ = *auxS; *auxS+= size_KZ*size_KZ;

  data->iXKZPivots = (PetscBLASInt*)*auxS;

  *auxS += FromIntToScalar(size_KZ);

  ierr = SlepcDenseCopy(data->iXKZ,data->ldiXKZ,data->XKZ,data->ldXKZ,size_KZ,size_KZ);CHKERRQ(ierr);

  ldXKZ = PetscBLASIntCast(data->ldiXKZ);

  LAPACKgetrf_(&s, &s, data->iXKZ, &ldXKZ, data->iXKZPivots, &info);

  if (info) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_LIB, "Error in Lapack XGETRF %d", info);

  PetscFunctionReturn(0);

}

#define DVD_COMPLEX_RAYLEIGH_QUOTIENT(ur,ui,Axr,Axi,Bxr,Bxi,eigr,eigi,b,ierr)\

{ \

  ierr = VecDot((Axr), (ur), &(b)[0]); CHKERRQ(ierr); /* r*A*r */ \

  ierr = VecDot((Axr), (ui), &(b)[1]); CHKERRQ(ierr); /* i*A*r */ \

  ierr = VecDot((Axi), (ur), &(b)[2]); CHKERRQ(ierr); /* r*A*i */ \

  ierr = VecDot((Axi), (ui), &(b)[3]); CHKERRQ(ierr); /* i*A*i */ \

  ierr = VecDot((Bxr), (ur), &(b)[4]); CHKERRQ(ierr); /* r*B*r */ \

  ierr = VecDot((Bxr), (ui), &(b)[5]); CHKERRQ(ierr); /* i*B*r */ \

  ierr = VecDot((Bxi), (ur), &(b)[6]); CHKERRQ(ierr); /* r*B*i */ \

  ierr = VecDot((Bxi), (ui), &(b)[7]); CHKERRQ(ierr); /* i*B*i */ \

  (b)[0]  = (b)[0]+(b)[3]; /* rAr+iAi */ \

  (b)[2] =  (b)[2]-(b)[1]; /* rAi-iAr */ \

  (b)[4] = (b)[4]+(b)[7]; /* rBr+iBi */ \

  (b)[6] = (b)[6]-(b)[5]; /* rBi-iBr */ \

  (b)[7] = (b)[4]*(b)[4] + (b)[6]*(b)[6]; /* k */ \

  *(eigr) = ((b)[0]*(b)[4] + (b)[2]*(b)[6]) / (b)[7]; /* eig_r */ \

  *(eigi) = ((b)[2]*(b)[4] - (b)[0]*(b)[6]) / (b)[7]; /* eig_i */ \

}