Subversion Repositories slepc-dev

Rev

Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1619 slepc 1 
/*
2 
  SLEPc eigensolver: "davidson"
3 
 
4 
  Step: calc the best eigenpairs in the subspace V.
5 
 
6 
  For that, performs these steps:
7 
    1) Update W <- A * V
8 
    2) Update H <- V' * W
9 
    3) Obtain eigenpairs of H
10 
    4) Select some eigenpairs
11 
    5) Compute the Ritz pairs of the selected ones
2110 jroman 12 
 
13 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
14 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2575 eromero 15 
   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
2110 jroman 16 
 
17 
   This file is part of SLEPc.
18 
 
19 
   SLEPc is free software: you can redistribute it and/or modify it under  the
20 
   terms of version 3 of the GNU Lesser General Public License as published by
21 
   the Free Software Foundation.
22 
 
23 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
24 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
25 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
26 
   more details.
27 
 
28 
   You  should have received a copy of the GNU Lesser General  Public  License
29 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
30 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1619 slepc 31 
*/
32 
 
1985 eromero 33 
#include "davidson.h"
2283 jroman 34 
#include <slepcblaslapack.h>
1619 slepc 35 
 
2555 eromero 36 
PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d);
2244 eromero 37 
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d);
2555 eromero 38 
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d);
2779 eromero 39 
PetscErrorCode dvd_calcpairs_projeig_ind(dvdDashboard *d);
2021 eromero 40 
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d);
41 
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d);
42 
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n);
2555 eromero 43 
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n);
2021 eromero 44 
PetscErrorCode dvd_calcpairs_X(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
45 
                               Vec *X);
46 
PetscErrorCode dvd_calcpairs_Y(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
47 
                               Vec *Y);
48 
PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
2605 eromero 49 
                                   Vec *R);
2682 jroman 50 
PetscErrorCode dvd_calcpairs_eig_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R);
2021 eromero 51 
PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d, PetscInt r_s,
52 
                                      PetscInt r_e, Vec *R);
2605 eromero 53 
PetscErrorCode dvd_calcpairs_updateV0(dvdDashboard *d, DvdReduction *r,
54 
                                      DvdMult_copy_func **sr);
55 
PetscErrorCode dvd_calcpairs_updateV1(dvdDashboard *d);
56 
PetscErrorCode dvd_calcpairs_updateW0(dvdDashboard *d, DvdReduction *r,
57 
                                      DvdMult_copy_func **sr);
58 
PetscErrorCode dvd_calcpairs_updateW1(dvdDashboard *d);
59 
PetscErrorCode dvd_calcpairs_updateAV0(dvdDashboard *d);
60 
PetscErrorCode dvd_calcpairs_updateAV1(dvdDashboard *d, DvdReduction *r,
61 
                                       DvdMult_copy_func **sr);
62 
PetscErrorCode dvd_calcpairs_updateBV0(dvdDashboard *d);
63 
PetscErrorCode dvd_calcpairs_updateBV1(dvdDashboard *d, DvdReduction *r,
64 
                                       DvdMult_copy_func **sr);
65 
PETSC_STATIC_INLINE PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,Vec *real_BV,PetscInt *size_cX,Vec **BV,PetscInt *size_BV,PetscInt *max_size_BV,PetscBool BV_shift,PetscInt *cX_in_proj,PetscScalar *MTX);
1619 slepc 66 
 
1735 eromero 67 
/**** Control routines ********************************************************/
1750 eromero 68 
#undef __FUNCT__  
69 
#define __FUNCT__ "dvd_calcpairs_qz"
2605 eromero 70 
PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d, dvdBlackboard *b,
71 
                                orthoV_type_t orth, IP ipI,
72 
                                PetscInt cX_proj, PetscBool harm)
1750 eromero 73 
{
2555 eromero 74 
  PetscErrorCode  ierr;
2654 jroman 75 
  PetscInt        i,max_cS;
2779 eromero 76 
  PetscBool       std_probl,her_probl,ind_probl;
1750 eromero 77 
 
78 
  PetscFunctionBegin;
79 
 
2012 eromero 80 
  std_probl = DVD_IS(d->sEP, DVD_EP_STD)?PETSC_TRUE:PETSC_FALSE;
2555 eromero 81 
  her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
2779 eromero 82 
  ind_probl = DVD_IS(d->sEP, DVD_EP_INDEFINITE)?PETSC_TRUE:PETSC_FALSE;
1750 eromero 83 
 
84 
  /* Setting configuration constrains */
2320 jroman 85 
#if !defined(PETSC_USE_COMPLEX)
1966 eromero 86 
  /* if the last converged eigenvalue is complex its conjugate pair is also
87 
     converged */
2605 eromero 88 
  b->max_nev = PetscMax(b->max_nev, d->nev+(her_probl && !d->B?0:1));
1966 eromero 89 
#else
90 
  b->max_nev = PetscMax(b->max_nev, d->nev);
91 
#endif
2605 eromero 92 
  b->max_size_proj = PetscMax(b->max_size_proj, b->max_size_V+cX_proj);
93 
  d->size_real_V = b->max_size_V+b->max_nev;
94 
  d->W_shift = d->B?PETSC_TRUE:PETSC_FALSE;
95 
  d->size_real_W = harm?(b->max_size_V+(d->W_shift?b->max_nev:b->max_size_cP)):0;
96 
  d->size_real_AV = b->max_size_V+b->max_size_cP;
97 
  d->size_BDS = 0;
2779 eromero 98 
  if (d->B && (her_probl || ind_probl) && (orth == DVD_ORTHOV_I || orth == DVD_ORTHOV_BOneMV)) {
2605 eromero 99 
    d->size_real_BV = b->size_V; d->BV_shift = PETSC_TRUE;
100 
    if (orth == DVD_ORTHOV_BOneMV) d->size_BDS = d->eps->nds;
101 
  } else if (d->B) {
102 
    d->size_real_BV = b->max_size_V + b->max_size_P; d->BV_shift = PETSC_FALSE;
103 
  } else {
104 
    d->size_real_BV = 0; d->BV_shift = PETSC_FALSE;
105 
  }
106 
  b->own_vecs+= d->size_real_V + d->size_real_W + d->size_real_AV +
107 
                d->size_real_BV + d->size_BDS;
108 
  b->own_scalars+= b->max_size_proj*b->max_size_proj*2*(std_probl?1:2) +
1757 eromero 109 
                                              /* H, G?, S, T? */
2605 eromero 110 
                   b->max_size_proj*b->max_size_proj*(std_probl?1:2) +
1757 eromero 111 
                                              /* pX, pY? */
2779 eromero 112 
                   b->max_nev*b->max_nev*(her_probl?0:(!d->B?1:2)) +
2605 eromero 113 
                                                /* cS?, cT? */
2779 eromero 114 
                   FromRealToScalar(d->size_real_V)*(ind_probl?1:0) + /* nBV */
115 
                   FromRealToScalar(b->max_size_proj)*(ind_probl?1:0); /* nBpX */
2605 eromero 116 
  b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X);
117 
                                                /* updateV0 */
2654 jroman 118 
  max_cS = PetscMax(b->max_size_X,cX_proj);
1820 eromero 119 
  b->max_size_auxS = PetscMax(PetscMax(
2605 eromero 120 
    b->max_size_auxS,
2654 jroman 121 
    b->max_size_X*b->max_size_proj*2*(std_probl?1:2) + /* updateAV1,BV1 */
122 
      max_cS*b->max_nev*(her_probl?0:(!d->B?1:2)) + /* updateV0,W0 */
2605 eromero 123 
                                                     /* SlepcReduction: in */
124 
      PetscMax(
2654 jroman 125 
        b->max_size_X*b->max_size_proj*2*(std_probl?1:2) + /* updateAV1,BV1 */
126 
          max_cS*b->max_nev*(her_probl?0:(!d->B?1:2)), /* updateV0,W0 */
2605 eromero 127 
                                                    /* SlepcReduction: out */
128 
        PetscMax(
129 
          b->max_size_proj*b->max_size_proj, /* updateAV0,BV0 */
130 
          b->max_size_proj+b->max_nev))), /* dvd_orth */
131 
    std_probl?0:(b->max_size_proj*11+16) /* projeig */);
2556 eromero 132 
#if defined(PETSC_USE_COMPLEX)
133 
  b->max_size_auxS = PetscMax(b->max_size_auxS, b->max_size_V);
134 
                                           /* dvd_calcpairs_projeig_eig */
135 
#endif
2779 eromero 136 
  if (ind_probl)
137 
    b->max_size_auxS = PetscMax(b->max_size_auxS,2*b->max_size_proj);
138 
                                           /* SlepcDenseIOrth */
1750 eromero 139 
 
140 
  /* Setup the step */
141 
  if (b->state >= DVD_STATE_CONF) {
2605 eromero 142 
    d->max_cX_in_proj = cX_proj;
143 
    d->max_size_P = b->max_size_P;
144 
    d->real_V = b->free_vecs; b->free_vecs+= d->size_real_V;
145 
    if (harm) {
146 
      d->real_W = b->free_vecs; b->free_vecs+= d->size_real_W;
147 
    } else {
148 
      d->real_W = PETSC_NULL;
149 
    }
150 
    d->real_AV = d->AV = b->free_vecs; b->free_vecs+= d->size_real_AV;
151 
    d->max_size_proj = b->max_size_proj;
152 
    d->real_H = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
153 
    d->ldH = b->max_size_proj;
154 
    d->pX = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
155 
    d->S = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
2555 eromero 156 
    if (!her_probl) {
157 
      d->cS = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
2605 eromero 158 
      d->max_size_cS = d->ldcS = b->max_nev;
2555 eromero 159 
    } else {
160 
      d->cS = PETSC_NULL;
2605 eromero 161 
      d->max_size_cS = d->ldcS = 0;
2779 eromero 162 
      d->orthoV_type = orth;
163 
      if (ind_probl) {
164 
        d->real_nBV = b->free_scalars; b->free_scalars+= FromRealToScalar(d->size_real_V);
165 
        d->nBpX = b->free_scalars; b->free_scalars+= FromRealToScalar(d->max_size_proj);
166 
      } else
167 
        d->real_nBV = d->nBDS = d->nBpX = PETSC_NULL;
2555 eromero 168 
    }
1795 eromero 169 
    d->ipV = ipI;
170 
    d->ipW = ipI;
2605 eromero 171 
    if (orth == DVD_ORTHOV_BOneMV) {
2555 eromero 172 
      d->BDS = b->free_vecs; b->free_vecs+= d->eps->nds;
173 
      for (i=0; i<d->eps->nds; i++) {
174 
        ierr = MatMult(d->B, d->eps->DS[i], d->BDS[i]); CHKERRQ(ierr);
175 
      }
2605 eromero 176 
    } else
177 
      d->BDS = PETSC_NULL;
1883 eromero 178 
    if (d->B) {
2605 eromero 179 
      d->real_BV = b->free_vecs; b->free_vecs+= d->size_real_BV;
2555 eromero 180 
    } else {
2605 eromero 181 
      d->size_real_BV = 0;
182 
      d->real_BV = PETSC_NULL;
183 
      d->BV_shift = PETSC_FALSE;
1883 eromero 184 
    }
185 
    if (!std_probl) {
2605 eromero 186 
      d->real_G = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
187 
      d->T = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
188 
      d->pY = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
189 
    } else {
190 
      d->real_G = PETSC_NULL;
191 
      d->T = PETSC_NULL;
192 
      d->pY = PETSC_NULL;
193 
    }
194 
    if (d->B && !her_probl) {
1966 eromero 195 
      d->cT = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
196 
      d->ldcT = b->max_nev;
2247 eromero 197 
    } else {
198 
      d->cT = PETSC_NULL;
199 
      d->ldcT = 0;
1750 eromero 200 
    }
201 
 
2555 eromero 202 
    d->calcPairs = dvd_calcpairs_proj;
1750 eromero 203 
    d->calcpairs_residual = dvd_calcpairs_res_0;
2682 jroman 204 
    d->calcpairs_residual_eig = dvd_calcpairs_eig_res_0;
1883 eromero 205 
    d->calcpairs_proj_res = dvd_calcpairs_proj_res;
2555 eromero 206 
    d->calcpairs_selectPairs = PETSC_NULL;
1751 eromero 207 
    d->ipI = ipI;
2244 eromero 208 
    DVD_FL_ADD(d->startList, dvd_calcpairs_qz_start);
1750 eromero 209 
  }
210 
 
211 
  PetscFunctionReturn(0);
212 
}
213 
 
214 
 
2244 eromero 215 
#undef __FUNCT__  
216 
#define __FUNCT__ "dvd_calcpairs_qz_start"
217 
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
218 
{
2779 eromero 219 
  PetscBool her_probl,ind_probl;
2654 jroman 220 
  PetscInt  i;
2244 eromero 221 
 
222 
  PetscFunctionBegin;
223 
  her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
2779 eromero 224 
  ind_probl = DVD_IS(d->sEP, DVD_EP_INDEFINITE)?PETSC_TRUE:PETSC_FALSE;
2244 eromero 225 
 
2605 eromero 226 
  d->size_V = 0;
227 
  d->V = d->real_V;
228 
  d->cX = d->real_V;
229 
  d->size_cX = 0;
230 
  d->max_size_V = d->size_real_V;
231 
  d->W = d->real_W;
232 
  d->max_size_W = d->size_real_W;
233 
  d->size_W = 0;
2244 eromero 234 
  d->size_AV = 0;
235 
  d->AV = d->real_AV;
2605 eromero 236 
  d->max_size_AV = d->size_real_AV;
2244 eromero 237 
  d->size_H = 0;
238 
  d->H = d->real_H;
2605 eromero 239 
  if (d->cS) for (i=0; i<d->max_size_cS*d->max_size_cS; i++) d->cS[i] = 0.0;
2244 eromero 240 
  d->size_BV = 0;
2605 eromero 241 
  d->BV = d->real_BV;
242 
  d->max_size_BV = d->size_real_BV;
2244 eromero 243 
  d->size_G = 0;
2247 eromero 244 
  d->G = d->real_G;
2605 eromero 245 
  if (d->cT) for (i=0; i<d->max_size_cS*d->max_size_cS; i++) d->cT[i] = 0.0;
246 
  d->cY = d->B && !her_probl ? d->W : PETSC_NULL;
247 
  d->BcX = d->orthoV_type == DVD_ORTHOV_I && d->B && her_probl ? d->BcX : PETSC_NULL;
248 
  d->size_cY = 0;
249 
  d->size_BcX = 0;
250 
  d->cX_in_V = d->cX_in_H = d->cX_in_G = d->cX_in_W = d->cX_in_AV = d->cX_in_BV = 0;
2779 eromero 251 
  d->nBV = d->nBcX = d->real_nBV;
2244 eromero 252 
  PetscFunctionReturn(0);
253 
}
254 
 
255 
 
1750 eromero 256 
#undef __FUNCT__
2555 eromero 257 
#define __FUNCT__ "dvd_calcpairs_proj"
258 
PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
1750 eromero 259 
{
1753 eromero 260 
  PetscErrorCode  ierr;
261 
  DvdReduction    r;
2663 jroman 262 
#define MAX_OPS 7
1753 eromero 263 
  DvdReductionChunk
2605 eromero 264 
                  ops[MAX_OPS];
1753 eromero 265 
  DvdMult_copy_func
2605 eromero 266 
                  sr[MAX_OPS], *sr0 = sr;
267 
  PetscInt        size_in, i;
268 
  PetscScalar     *in = d->auxS, *out;
269 
  PetscBool       stdp;
1753 eromero 270 
 
1750 eromero 271 
  PetscFunctionBegin;
272 
 
2605 eromero 273 
  stdp = DVD_IS(d->sEP, DVD_EP_STD)?PETSC_TRUE:PETSC_FALSE;
274 
  size_in =
2654 jroman 275 
    (d->size_cX+d->V_tra_s-d->cX_in_H)*d->V_tra_s*(d->cT?2:(d->cS?1:0)) + /* updateV0,W0 */
2605 eromero 276 
    (d->size_H*(d->V_new_e-d->V_new_s)*2+
2654 jroman 277 
      (d->V_new_e-d->V_new_s)*(d->V_new_e-d->V_new_s))*(!stdp?2:1); /* updateAV1,BV1 */
2605 eromero 278 
 
279 
  out = in+size_in;
280 
 
2682 jroman 281 
  /* Check consistency */
282 
  if (2*size_in > d->size_auxS) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
2605 eromero 283 
 
1753 eromero 284 
  /* Prepare reductions */
2605 eromero 285 
  ierr = SlepcAllReduceSumBegin(ops, MAX_OPS, in, out, size_in, &r,
1753 eromero 286 
                                ((PetscObject)d->V[0])->comm); CHKERRQ(ierr);
2605 eromero 287 
  /* Allocate size_in */
288 
  d->auxS+= size_in;
289 
  d->size_auxS-= size_in;
1753 eromero 290 
 
2555 eromero 291 
  /* Update AV, BV, W and the projected matrices */
2605 eromero 292 
  /* 1. S <- S*MT */
293 
  ierr = dvd_calcpairs_updateV0(d, &r, &sr0); CHKERRQ(ierr);
294 
  ierr = dvd_calcpairs_updateW0(d, &r, &sr0); CHKERRQ(ierr);
295 
  ierr = dvd_calcpairs_updateAV0(d); CHKERRQ(ierr);
296 
  ierr = dvd_calcpairs_updateBV0(d); CHKERRQ(ierr);
297 
  /* 2. V <- orth(V, V_new) */
298 
  ierr = dvd_calcpairs_updateV1(d); CHKERRQ(ierr);
299 
  /* 3. AV <- [AV A * V(V_new_s:V_new_e-1)] */
300 
  /* Check consistency */
301 
  if (d->size_AV != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
302 
  for (i=d->V_new_s; i<d->V_new_e; i++) {
303 
    ierr = MatMult(d->A, d->V[i], d->AV[i]); CHKERRQ(ierr);
2555 eromero 304 
  }
2605 eromero 305 
  d->size_AV = d->V_new_e;
306 
  /* 4. BV <- [BV B * V(V_new_s:V_new_e-1)] */
307 
  if (d->B && d->orthoV_type != DVD_ORTHOV_BOneMV) {
308 
    /* Check consistency */
309 
    if (d->size_BV != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
310 
    for (i=d->V_new_s; i<d->V_new_e; i++) {
311 
      ierr = MatMult(d->B, d->V[i], d->BV[i]); CHKERRQ(ierr);
312 
    }
313 
    d->size_BV = d->V_new_e;
2021 eromero 314 
  }
2605 eromero 315 
  /* 5 <- W <- [W f(AV,BV)] */
316 
  ierr = dvd_calcpairs_updateW1(d); CHKERRQ(ierr);
317 
  ierr = dvd_calcpairs_updateAV1(d, &r, &sr0); CHKERRQ(ierr);
318 
  ierr = dvd_calcpairs_updateBV1(d, &r, &sr0); CHKERRQ(ierr);
1753 eromero 319 
 
2605 eromero 320 
  /* Deallocate size_in */
321 
  d->auxS-= size_in;
322 
  d->size_auxS+= size_in;
323 
 
1753 eromero 324 
  /* Do reductions */
325 
  ierr = SlepcAllReduceSumEnd(&r); CHKERRQ(ierr);
326 
 
2555 eromero 327 
  /* Perform the transformation on the projected problem */
328 
  if (d->calcpairs_proj_trans) {
329 
    ierr = d->calcpairs_proj_trans(d); CHKERRQ(ierr);
330 
  }
331 
 
2605 eromero 332 
  d->V_tra_s = d->V_tra_e = 0;
333 
  d->V_new_s = d->V_new_e;
2555 eromero 334 
 
335 
  /* Solve the projected problem */
2605 eromero 336 
  if (DVD_IS(d->sEP, DVD_EP_STD)) {
2779 eromero 337 
    if (DVD_IS(d->sEP, DVD_EP_INDEFINITE)) {
338 
      ierr = dvd_calcpairs_projeig_ind(d); CHKERRQ(ierr);
339 
    } else if (DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
2605 eromero 340 
      ierr = dvd_calcpairs_projeig_eig(d); CHKERRQ(ierr);
2021 eromero 341 
    } else {
2605 eromero 342 
      ierr = dvd_calcpairs_projeig_qz_std(d); CHKERRQ(ierr);
2021 eromero 343 
    }
2605 eromero 344 
  } else {
345 
    ierr = dvd_calcpairs_projeig_qz_gen(d); CHKERRQ(ierr);
346 
  }
1619 slepc 347 
 
1829 eromero 348 
  /* Check consistency */
2605 eromero 349 
  if (d->size_V != d->V_new_e || d->size_V+d->cX_in_H != d->size_H || d->cX_in_V != d->cX_in_H ||
350 
      d->size_V != d->size_AV || d->cX_in_H != d->cX_in_AV ||
351 
        (DVD_ISNOT(d->sEP, DVD_EP_STD) && (
352 
          d->size_V+d->cX_in_G != d->size_G || d->cX_in_H != d->cX_in_G ||
353 
          d->size_H != d->size_G || (d->BV && (
354 
            d->size_V != d->size_BV || d->cX_in_H != d->cX_in_BV)))) ||
355 
      (d->W && d->size_W != d->size_V)) {
2214 jroman 356 
    SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
1829 eromero 357 
  }
358 
 
1619 slepc 359 
  PetscFunctionReturn(0);
2663 jroman 360 
#undef MAX_OPS
1619 slepc 361 
}
362 
 
1735 eromero 363 
/**** Basic routines **********************************************************/
364 
 
365 
#undef __FUNCT__
2605 eromero 366 
#define __FUNCT__ "dvd_calcpairs_updateV0"
367 
/* auxV: V_tra_s, DvdMult_copy_func: 1 */
368 
PetscErrorCode dvd_calcpairs_updateV0(dvdDashboard *d, DvdReduction *r,
369 
                                      DvdMult_copy_func **sr)
1795 eromero 370 
{
2021 eromero 371 
  PetscErrorCode  ierr;
2779 eromero 372 
  PetscInt        rm,i;
2605 eromero 373 
 
374 
  PetscFunctionBegin;
375 
 
376 
  if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
377 
 
2779 eromero 378 
  /* Update nBcX and nBV */
379 
  if (d->nBcX && d->nBpX && d->nBV) {
380 
    d->nBV+= d->V_tra_s;
381 
    for (i=0; i<d->V_tra_s; i++) d->nBcX[d->size_cX+i] = d->nBpX[i];
382 
    for (i=d->V_tra_s; i<d->V_tra_e; i++) d->nBV[i-d->V_tra_s] = d->nBpX[i];
383 
  }
384 
 
2605 eromero 385 
  /* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
386 
  ierr = dvd_calcpairs_updateBV0_gen(d,d->real_V,&d->size_cX,&d->V,&d->size_V,&d->max_size_V,PETSC_TRUE,&d->cX_in_V,d->MTX);CHKERRQ(ierr);
387 
 
388 
  /* Udpate cS for standard problems */
2624 eromero 389 
  if (d->cS && !d->cT && !d->cY && (d->V_tra_s > d->max_cX_in_proj || d->size_cX >= d->nev)) {
2605 eromero 390 
    /* Check consistency */
391 
    if (d->size_cS+d->V_tra_s != d->size_cX) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
392 
 
393 
    /* auxV <- AV * MTX(0:V_tra_e-1) */
2654 jroman 394 
    rm = d->size_cX>=d->nev?0:d->max_cX_in_proj;
395 
    ierr = SlepcUpdateVectorsZ(d->auxV, 0.0, 1.0, d->AV-d->cX_in_AV, d->size_AV+d->cX_in_AV, d->MTX, d->ldMTX, d->size_MT, d->V_tra_s-rm); CHKERRQ(ierr);
2605 eromero 396 
 
397 
    /* cS(:, size_cS:) <- cX' * auxV */
398 
    ierr = VecsMultS(&d->cS[d->ldcS*d->size_cS], 0, d->ldcS, d->cX, 0, d->size_cX-rm, d->auxV, 0, d->V_tra_s-rm, r, (*sr)++); CHKERRQ(ierr);
399 
    d->size_cS+= d->V_tra_s-rm;
400 
  }
2779 eromero 401 
 
2605 eromero 402 
  PetscFunctionReturn(0);
403 
}
404 
 
405 
 
406 
#undef __FUNCT__
407 
#define __FUNCT__ "dvd_calcpairs_updateV1"
408 
/* auxS: size_cX+V_new_e+1 */
409 
PetscErrorCode dvd_calcpairs_updateV1(dvdDashboard *d)
410 
{
411 
  PetscErrorCode  ierr;
1829 eromero 412 
  Vec             *cX = d->BcX? d->BcX : ( (d->cY && !d->W)? d->cY : d->cX );
1795 eromero 413 
 
414 
  PetscFunctionBegin;
415 
 
2605 eromero 416 
  if (d->V_new_s == d->V_new_e) { PetscFunctionReturn(0); }
417 
 
418 
  /* Check consistency */
419 
  if (d->size_V != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
420 
 
421 
  /* V <- gs([cX V(0:V_new_s-1)], V(V_new_s:V_new_e-1)) */
422 
  if (d->orthoV_type == DVD_ORTHOV_BOneMV) {
2779 eromero 423 
    ierr = dvd_BorthV(d->ipV, d->eps->DS, d->BDS, d->nBDS, d->eps->nds, d->cX, d->real_BV, d->nBcX,
424 
                      d->size_cX, d->V, d->BV, d->nBV, d->V_new_s, d->V_new_e,
2605 eromero 425 
                      d->auxS, d->eps->rand); CHKERRQ(ierr);
426 
    d->size_BV = d->V_new_e;
2555 eromero 427 
  } else {
428 
    ierr = dvd_orthV(d->ipV, d->eps->DS, d->eps->nds, cX, d->size_cX, d->V,
2605 eromero 429 
                   d->V_new_s, d->V_new_e, d->auxS, d->eps->rand);
2555 eromero 430 
    CHKERRQ(ierr);
431 
  }
2605 eromero 432 
  d->size_V = d->V_new_e;
1795 eromero 433 
 
2021 eromero 434 
  PetscFunctionReturn(0);
1795 eromero 435 
}
436 
 
437 
#undef __FUNCT__
2605 eromero 438 
#define __FUNCT__ "dvd_calcpairs_updateW0"
439 
/* auxV: V_tra_s, DvdMult_copy_func: 2 */
440 
PetscErrorCode dvd_calcpairs_updateW0(dvdDashboard *d, DvdReduction *r,
441 
                                      DvdMult_copy_func **sr)
1795 eromero 442 
{
2021 eromero 443 
  PetscErrorCode  ierr;
2605 eromero 444 
  PetscInt        rm;
1795 eromero 445 
 
446 
  PetscFunctionBegin;
447 
 
2682 jroman 448 
  if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
2605 eromero 449 
 
450 
  /* cY <- [cY W*MTY(0:V_tra_s-1)], W <- W*MTY(V_tra_s:V_tra_e) */
451 
  ierr = dvd_calcpairs_updateBV0_gen(d,d->real_W,&d->size_cY,&d->W,&d->size_W,&d->max_size_W,d->W_shift,&d->cX_in_W,d->MTY);CHKERRQ(ierr);
452 
 
453 
  /* Udpate cS and cT */
2682 jroman 454 
  if (d->cT && (d->V_tra_s > d->max_cX_in_proj || d->size_cX >= d->nev)) {
2605 eromero 455 
    /* Check consistency */
2682 jroman 456 
    if (d->size_cS+d->V_tra_s != d->size_cX || (d->W && d->size_cY != d->size_cX)) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
2605 eromero 457 
 
458 
    /* auxV <- AV * MTX(0:V_tra_e-1) */
459 
    rm = d->size_cX>=d->nev?0:d->max_cX_in_proj;
2654 jroman 460 
    ierr = SlepcUpdateVectorsZ(d->auxV, 0.0, 1.0, d->AV-d->cX_in_H, d->size_AV-d->cX_in_H, d->MTX, d->ldMTX, d->size_MT, d->V_tra_s-rm); CHKERRQ(ierr);
2605 eromero 461 
 
462 
    /* cS(:, size_cS:) <- cY' * auxV */
2682 jroman 463 
    ierr = VecsMultS(&d->cS[d->ldcS*d->size_cS], 0, d->ldcS, d->cY?d->cY:d->cX, 0, d->size_cX-rm, d->auxV, 0, d->V_tra_s-rm, r, (*sr)++); CHKERRQ(ierr);
2605 eromero 464 
 
465 
    /* auxV <- BV * MTX(0:V_tra_e-1) */
2654 jroman 466 
    ierr = SlepcUpdateVectorsZ(d->auxV, 0.0, 1.0, d->BV-d->cX_in_H, d->size_BV-d->cX_in_H, d->MTX, d->ldMTX, d->size_MT, d->V_tra_s-rm); CHKERRQ(ierr);
2605 eromero 467 
 
468 
    /* cT(:, size_cS:) <- cY' * auxV */
2682 jroman 469 
    ierr = VecsMultS(&d->cT[d->ldcS*d->size_cS], 0, d->ldcS, d->cY?d->cY:d->cX, 0, d->size_cX-rm, d->auxV, 0, d->V_tra_s-rm, r, (*sr)++); CHKERRQ(ierr);
2654 jroman 470 
 
2605 eromero 471 
    d->size_cS+= d->V_tra_s-rm;
472 
    d->size_cT+= d->V_tra_s-rm;
473 
  }
474 
 
475 
  PetscFunctionReturn(0);
476 
}
477 
 
478 
 
479 
#undef __FUNCT__
480 
#define __FUNCT__ "dvd_calcpairs_updateW1"
481 
/* auxS: size_cX+V_new_e+1 */
482 
PetscErrorCode dvd_calcpairs_updateW1(dvdDashboard *d)
483 
{
484 
  PetscErrorCode  ierr;
2664 jroman 485 
  Vec             *cY = d->cY?d->cY:d->cX;
2605 eromero 486 
 
487 
  PetscFunctionBegin;
488 
 
2664 jroman 489 
  if (!d->W || d->V_new_s == d->V_new_e) { PetscFunctionReturn(0); }
2605 eromero 490 
 
491 
  /* Check consistency */
492 
  if (d->size_W != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
493 
 
1795 eromero 494 
  /* Update W */
2021 eromero 495 
  ierr = d->calcpairs_W(d); CHKERRQ(ierr);
1795 eromero 496 
 
2605 eromero 497 
  /* W <- gs([cY W(0:V_new_s-1)], W(V_new_s:V_new_e-1)) */
2664 jroman 498 
  ierr = dvd_orthV(d->ipW, PETSC_NULL, 0, cY, d->size_cX, d->W-d->cX_in_W, d->V_new_s+d->cX_in_W, d->V_new_e+d->cX_in_W, d->auxS, d->eps->rand); CHKERRQ(ierr);
2605 eromero 499 
  d->size_W = d->V_new_e;
1795 eromero 500 
 
2021 eromero 501 
  PetscFunctionReturn(0);
1795 eromero 502 
}
503 
 
504 
#undef __FUNCT__
2605 eromero 505 
#define __FUNCT__ "dvd_calcpairs_updateAV0"
506 
/* auxS: size_H*(V_tra_e-V_tra_s) */
507 
PetscErrorCode dvd_calcpairs_updateAV0(dvdDashboard *d)
1619 slepc 508 
{
2021 eromero 509 
  PetscErrorCode  ierr;
2605 eromero 510 
  PetscScalar     *MTY = d->W?d->MTY:d->MTX;
2654 jroman 511 
  PetscInt        cMT, tra_s;
1619 slepc 512 
 
513 
  PetscFunctionBegin;
514 
 
2605 eromero 515 
  if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
1619 slepc 516 
 
2605 eromero 517 
  /* AV(V_tra_s-cp-1:) = cAV*MTX(V_tra_s:) */
518 
  ierr = dvd_calcpairs_updateBV0_gen(d,d->real_AV,PETSC_NULL,&d->AV,&d->size_AV,&d->max_size_AV,PETSC_FALSE,&d->cX_in_AV,d->MTX);CHKERRQ(ierr);
519 
  tra_s = PetscMax(d->V_tra_s-d->max_cX_in_proj, 0);
520 
  cMT = d->V_tra_e - tra_s;
521 
 
522 
  /* Update H <- MTY(tra_s)' * (H * MTX(tra_s:)) */
523 
  ierr = SlepcDenseMatProdTriang(d->auxS, 0, d->ldH, d->H, d->sH, d->ldH, d->size_H, d->size_H, PETSC_FALSE, &d->MTX[d->ldMTX*tra_s], 0, d->ldMTX, d->size_MT, cMT, PETSC_FALSE); CHKERRQ(ierr);
524 
  ierr = SlepcDenseMatProdTriang(d->H, d->sH, d->ldH, &MTY[d->ldMTX*tra_s], 0, d->ldMTX, d->size_MT, cMT, PETSC_TRUE, d->auxS, 0, d->ldH, d->size_H, cMT, PETSC_FALSE); CHKERRQ(ierr);
525 
  d->size_H = cMT;
2654 jroman 526 
  d->cX_in_H = d->cX_in_AV;
2605 eromero 527 
 
2021 eromero 528 
  PetscFunctionReturn(0);
1735 eromero 529 
}
1640 slepc 530 
 
2605 eromero 531 
 
1735 eromero 532 
#undef __FUNCT__
2605 eromero 533 
#define __FUNCT__ "dvd_calcpairs_updateAV1"
534 
/* DvdMult_copy_func: 2 */
535 
PetscErrorCode dvd_calcpairs_updateAV1(dvdDashboard *d, DvdReduction *r,
536 
                                       DvdMult_copy_func **sr)
1735 eromero 537 
{
2021 eromero 538 
  PetscErrorCode  ierr;
2605 eromero 539 
  Vec             *W = d->W?d->W:d->V;
1619 slepc 540 
 
1735 eromero 541 
  PetscFunctionBegin;
1619 slepc 542 
 
2605 eromero 543 
  if (d->V_new_s == d->V_new_e) { PetscFunctionReturn(0); }
1619 slepc 544 
 
2605 eromero 545 
  /* Check consistency */
546 
  if (d->size_H != d->V_new_s+d->cX_in_H || d->size_V != d->V_new_e) {
547 
    SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
548 
  }
549 
 
550 
  /* H = [H               W(old)'*AV(new);
551 
          W(new)'*AV(old) W(new)'*AV(new) ],
552 
     being old=0:V_new_s-1, new=V_new_s:V_new_e-1 */
2624 eromero 553 
  ierr = VecsMultS(d->H,d->sH,d->ldH,W-d->cX_in_H,d->V_new_s+d->cX_in_H, d->V_new_e+d->cX_in_H, d->AV-d->cX_in_H,d->V_new_s+d->cX_in_H,d->V_new_e+d->cX_in_H, r, (*sr)++); CHKERRQ(ierr);
2605 eromero 554 
  d->size_H = d->V_new_e+d->cX_in_H;
555 
 
2021 eromero 556 
  PetscFunctionReturn(0);
1630 slepc 557 
}
1619 slepc 558 
 
1630 slepc 559 
#undef __FUNCT__
2605 eromero 560 
#define __FUNCT__ "dvd_calcpairs_updateBV0"
561 
/* auxS: max(BcX*(size_cX+V_new_e+1), size_G*(V_tra_e-V_tra_s)) */
562 
PetscErrorCode dvd_calcpairs_updateBV0(dvdDashboard *d)
1630 slepc 563 
{
2021 eromero 564 
  PetscErrorCode  ierr;
2605 eromero 565 
  PetscScalar     *MTY = d->W?d->MTY:d->MTX;
2654 jroman 566 
  PetscInt        cMT, tra_s, i;
2605 eromero 567 
  PetscBool       lindep;
568 
  PetscReal       norm;
1626 slepc 569 
 
1630 slepc 570 
  PetscFunctionBegin;
1619 slepc 571 
 
2605 eromero 572 
  if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
1619 slepc 573 
 
2605 eromero 574 
  /* BV <- BV*MTX */
575 
  ierr = dvd_calcpairs_updateBV0_gen(d,d->real_BV,PETSC_NULL,&d->BV,&d->size_BV,&d->max_size_BV,d->BV_shift,&d->cX_in_BV,d->MTX);CHKERRQ(ierr);
576 
 
577 
  /* If BcX, BcX <- orth(BcX) */
578 
  if (d->BcX) {
579 
    for (i=0; i<d->V_tra_s; i++) {
580 
      ierr = IPOrthogonalize(d->ipI, 0, PETSC_NULL, d->size_BcX+i, PETSC_NULL,
581 
                             d->BcX, d->BcX[d->size_BcX+i], PETSC_NULL,
582 
                             &norm, &lindep); CHKERRQ(ierr);
583 
      if(lindep) SETERRQ(((PetscObject)d->ipI)->comm,1, "Error during orth(BcX, B*cX(new))");
584 
      ierr = VecScale(d->BcX[d->size_BcX+i], 1.0/norm); CHKERRQ(ierr);
585 
    }
586 
    d->size_BcX+= d->V_tra_s;
587 
  }
588 
 
589 
  /* Update G <- MTY' * (G * MTX) */
590 
  if (d->G) {
591 
    tra_s = PetscMax(d->V_tra_s-d->max_cX_in_proj, 0);
592 
    cMT = d->V_tra_e - tra_s;
593 
    ierr = SlepcDenseMatProdTriang(d->auxS, 0, d->ldH, d->G, d->sG, d->ldH, d->size_G, d->size_G, PETSC_FALSE, &d->MTX[d->ldMTX*tra_s], 0, d->ldMTX, d->size_MT, cMT, PETSC_FALSE); CHKERRQ(ierr);
594 
    ierr = SlepcDenseMatProdTriang(d->G, d->sG, d->ldH, &MTY[d->ldMTX*tra_s], 0, d->ldMTX, d->size_MT, cMT, PETSC_TRUE, d->auxS, 0, d->ldH, d->size_G, cMT, PETSC_FALSE); CHKERRQ(ierr);
595 
    d->size_G = cMT;
2654 jroman 596 
    d->cX_in_G = d->cX_in_V;
2605 eromero 597 
  }
598 
 
2021 eromero 599 
  PetscFunctionReturn(0);
1619 slepc 600 
}
601 
 
1795 eromero 602 
 
1630 slepc 603 
#undef __FUNCT__
2605 eromero 604 
#define __FUNCT__ "dvd_calcpairs_updateBV1"
605 
/* DvdMult_copy_func: 2 */
606 
PetscErrorCode dvd_calcpairs_updateBV1(dvdDashboard *d, DvdReduction *r,
607 
                                       DvdMult_copy_func **sr)
1619 slepc 608 
{
2021 eromero 609 
  PetscErrorCode  ierr;
2605 eromero 610 
  Vec             *W = d->W?d->W:d->V, *BV = d->BV?d->BV:d->V;
1619 slepc 611 
 
612 
  PetscFunctionBegin;
613 
 
2605 eromero 614 
  if (!d->G || d->V_new_s == d->V_new_e) { PetscFunctionReturn(0); }
1630 slepc 615 
 
2605 eromero 616 
  /* G = [G               W(old)'*BV(new);
617 
          W(new)'*BV(old) W(new)'*BV(new) ],
618 
     being old=0:V_new_s-1, new=V_new_s:V_new_e-1 */
2624 eromero 619 
  ierr = VecsMultS(d->G,d->sG,d->ldH,W-d->cX_in_G,d->V_new_s+d->cX_in_G,d->V_new_e+d->cX_in_G,BV-d->cX_in_G,d->V_new_s+d->cX_in_G,d->V_new_e+d->cX_in_G,r,(*sr)++); CHKERRQ(ierr);
2605 eromero 620 
  d->size_G = d->V_new_e+d->cX_in_G;
621 
 
2021 eromero 622 
  PetscFunctionReturn(0);
1619 slepc 623 
}
624 
 
2556 eromero 625 
/* in complex, d->size_H real auxiliar values are needed */
1735 eromero 626 
#undef __FUNCT__
2555 eromero 627 
#define __FUNCT__ "dvd_calcpairs_projeig_eig"
628 
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d)
629 
{
630 
  PetscErrorCode  ierr;
2556 eromero 631 
  PetscReal       *w;
632 
#if defined(PETSC_USE_COMPLEX)
633 
  PetscInt        i;
634 
#endif
2555 eromero 635 
 
636 
  PetscFunctionBegin;
637 
 
638 
  /* S <- H */
639 
  d->ldS = d->ldpX = d->size_H;
640 
  ierr = SlepcDenseCopyTriang(d->S, DVD_MAT_LTRIANG, d->size_H, d->H, d->sH, d->ldH,
2605 eromero 641 
                              d->size_H, d->size_H); CHKERRQ(ierr);
2555 eromero 642 
 
643 
  /* S = pX' * L * pX */
2556 eromero 644 
#if !defined(PETSC_USE_COMPLEX)
2605 eromero 645 
  w = d->eigr-d->cX_in_H;
2556 eromero 646 
#else
647 
  w = (PetscReal*)d->auxS;
648 
#endif
649 
  ierr = EPSDenseHEP(d->size_H, d->S, d->ldS, w, d->pX); CHKERRQ(ierr);
650 
#if defined(PETSC_USE_COMPLEX)
2605 eromero 651 
  for (i=0; i<d->size_H; i++) d->eigr[i-d->cX_in_H] = w[i];
2556 eromero 652 
#endif
2555 eromero 653 
 
654 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
655 
 
656 
  PetscFunctionReturn(0);
657 
}
658 
 
659 
 
2779 eromero 660 
/* in complex, d->size_H real auxiliar values are needed */
2555 eromero 661 
#undef __FUNCT__
2779 eromero 662 
#define __FUNCT__ "dvd_calcpairs_projeig_ind"
663 
PetscErrorCode dvd_calcpairs_projeig_ind(dvdDashboard *d)
664 
{
665 
  PetscErrorCode  ierr;
666 
  PetscInt        i,j;
667 
 
668 
  PetscFunctionBegin;
669 
 
670 
  /* S <- H */
671 
  d->ldS = d->ldpX = d->size_H;
672 
  ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
673 
                              d->size_H, d->size_H); CHKERRQ(ierr);
674 
 
675 
  /* S <- nBV * S */
676 
  for (i=0; i<d->size_H; i++) {
677 
    if (d->nBV[i] < 0.0)
678 
      for (j=0; j<d->size_H; j++) d->S[i+j*d->ldS]*= d->nBV[i];
679 
  }
680 
 
681 
  /* S = pX * L */
682 
  ierr = EPSDenseNHEP(d->size_H, d->S, d->eigr-d->cX_in_H, d->eigi-d->cX_in_H, d->pX, PETSC_NULL); CHKERRQ(ierr);
683 
 
684 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
685 
 
686 
  PetscFunctionReturn(0);
687 
}
688 
 
689 
#undef __FUNCT__
1750 eromero 690 
#define __FUNCT__ "dvd_calcpairs_projeig_qz_std"
2021 eromero 691 
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d)
1636 slepc 692 
{
693 
  PetscErrorCode  ierr;
694 
 
695 
  PetscFunctionBegin;
696 
 
1750 eromero 697 
  /* S <- H */
698 
  d->ldS = d->ldpX = d->size_H;
1752 eromero 699 
  ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
700 
                              d->size_H, d->size_H);
1750 eromero 701 
 
702 
  /* S = pX' * H * pX */
2021 eromero 703 
  ierr = EPSDenseHessenberg(d->size_H, 0, d->S, d->ldS, d->pX); CHKERRQ(ierr);
2605 eromero 704 
  ierr = EPSDenseSchur(d->size_H, 0, d->S, d->ldS, d->pX, d->eigr-d->cX_in_H, d->eigi-d->cX_in_H);
1750 eromero 705 
  CHKERRQ(ierr);
1636 slepc 706 
 
2555 eromero 707 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_qz;
708 
 
1636 slepc 709 
  PetscFunctionReturn(0);
710 
}
711 
 
2223 jroman 712 
#undef __FUNCT__
713 
#define __FUNCT__ "dvd_calcpairs_projeig_qz_gen"
1750 eromero 714 
/*
715 
  auxS(dgges) = size_H (beta) + 8*size_H+16 (work)
716 
  auxS(zgges) = size_H (beta) + 1+2*size_H (work) + 8*size_H (rwork)
717 
*/
2021 eromero 718 
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d)
1750 eromero 719 
{
2392 jroman 720 
#if defined(SLEPC_MISSING_LAPACK_GGES)
721 
  PetscFunctionBegin;
722 
  SETERRQ(((PetscObject)(d->eps))->comm,PETSC_ERR_SUP,"GGES - Lapack routine is unavailable.");
723 
#else
1750 eromero 724 
  PetscErrorCode  ierr;
725 
  PetscScalar     *beta = d->auxS;
726 
#if !defined(PETSC_USE_COMPLEX)
727 
  PetscScalar     *auxS = beta + d->size_H;
728 
  PetscBLASInt    n_auxS = d->size_auxS - d->size_H;
729 
#else
1883 eromero 730 
  PetscReal       *auxR = (PetscReal*)(beta + d->size_H);
1750 eromero 731 
  PetscScalar     *auxS = (PetscScalar*)(auxR+8*d->size_H);
1883 eromero 732 
  PetscBLASInt    n_auxS = d->size_auxS - 9*d->size_H;
1750 eromero 733 
#endif
1795 eromero 734 
  PetscInt        i;
1750 eromero 735 
  PetscBLASInt    info,n,a;
736 
 
737 
  PetscFunctionBegin;
2654 jroman 738 
  if (d->pY) { PetscValidScalarPointer(d->pY,0); }
739 
  PetscValidScalarPointer(d->S,0);
740 
  PetscValidScalarPointer(d->T,0);
741 
  PetscValidScalarPointer(d->eigr,0);
742 
  PetscValidScalarPointer(d->eigi,0);
743 
  PetscValidScalarPointer(d->auxS,0);
744 
 
1750 eromero 745 
  /* S <- H, T <- G */
746 
  d->ldS = d->ldT = d->ldpX = d->ldpY = d->size_H;
1752 eromero 747 
  ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
2139 jroman 748 
                              d->size_H, d->size_H);CHKERRQ(ierr);
1752 eromero 749 
  ierr = SlepcDenseCopyTriang(d->T, 0, d->size_H, d->G, d->sG, d->ldH,
2139 jroman 750 
                              d->size_H, d->size_H);CHKERRQ(ierr);
1750 eromero 751 
 
752 
  /* S = Z'*H*Q, T = Z'*G*Q */
753 
  n = d->size_H;
754 
#if !defined(PETSC_USE_COMPLEX)
1757 eromero 755 
  LAPACKgges_(d->pY?"V":"N", "V", "N", PETSC_NULL, &n, d->S, &n, d->T, &n,
2624 eromero 756 
              &a, d->eigr-d->cX_in_H, d->eigi-d->cX_in_H, beta, d->pY, &n, d->pX, &n,
1750 eromero 757 
              auxS, &n_auxS, PETSC_NULL, &info);
758 
#else
1757 eromero 759 
  LAPACKgges_(d->pY?"V":"N", "V", "N", PETSC_NULL, &n, d->S, &n, d->T, &n,
2624 eromero 760 
              &a, d->eigr-d->cX_in_H, beta, d->pY, &n, d->pX, &n,
1750 eromero 761 
              auxS, &n_auxS, auxR, PETSC_NULL, &info);
762 
#endif
2214 jroman 763 
  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB, "Error in Lapack GGES %d", info);
1750 eromero 764 
 
765 
  /* eigr[i] <- eigr[i] / beta[i] */
2624 eromero 766 
  for (i=0; i<d->size_H; i++)
767 
    d->eigr[i-d->cX_in_H] /= beta[i],
768 
    d->eigi[i-d->cX_in_H] /= beta[i];
1750 eromero 769 
 
2555 eromero 770 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_qz;
771 
 
1750 eromero 772 
  PetscFunctionReturn(0);
2392 jroman 773 
#endif
1750 eromero 774 
}
775 
 
2555 eromero 776 
#undef __FUNCT__
777 
#define __FUNCT__ "dvd_calcpairs_selectPairs_eig"
778 
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n)
779 
{
780 
  PetscErrorCode  ierr;
1750 eromero 781 
 
2555 eromero 782 
  PetscFunctionBegin;
783 
 
2605 eromero 784 
  ierr = EPSSortDenseHEP(d->eps, d->size_H, 0, d->eigr-d->cX_in_H, d->pX, d->ldpX);
2555 eromero 785 
  CHKERRQ(ierr);
786 
 
787 
  if (d->calcpairs_eigs_trans) {
788 
    ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
789 
  }
790 
 
2779 eromero 791 
  if (DVD_IS(d->sEP, DVD_EP_INDEFINITE)) {
792 
    ierr = SlepcDenseIOrth(d->pX,d->ldpX,d->size_H,d->size_H,d->nBV-d->cX_in_H,d->nBpX,d->auxS,d->size_auxS);CHKERRQ(ierr);
793 
  }
794 
 
2555 eromero 795 
  PetscFunctionReturn(0);
796 
}
797 
 
798 
 
1750 eromero 799 
#undef __FUNCT__
800 
#define __FUNCT__ "dvd_calcpairs_selectPairs_qz"
2021 eromero 801 
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n)
1750 eromero 802 
{
1805 eromero 803 
  PetscErrorCode  ierr;
1750 eromero 804 
  PetscFunctionBegin;
805 
 
1805 eromero 806 
  if ((d->ldpX != d->size_H) ||
807 
      ( d->T &&
808 
        ((d->ldS != d->ldT) || (d->ldpX != d->ldpY) ||
809 
         (d->ldpX != d->size_H)) ) ) {
2214 jroman 810 
     SETERRQ(PETSC_COMM_SELF,1, "Error before ordering eigenpairs");
1750 eromero 811 
  }
812 
 
1829 eromero 813 
  if (d->T) {
2624 eromero 814 
    ierr = EPSSortDenseSchurGeneralized(d->eps, d->size_H, 0, n, d->S, d->T,
2605 eromero 815 
                                        d->ldS, d->pY, d->pX, d->eigr-d->cX_in_H,
816 
                                        d->eigi-d->cX_in_H); CHKERRQ(ierr);
1829 eromero 817 
  } else {
2624 eromero 818 
    ierr = EPSSortDenseSchur(d->eps, d->size_H, 0, d->S, d->ldS, d->pX,
2605 eromero 819 
                             d->eigr-d->cX_in_H, d->eigi-d->cX_in_H); CHKERRQ(ierr);
1829 eromero 820 
  }
1750 eromero 821 
 
2021 eromero 822 
  if (d->calcpairs_eigs_trans) {
823 
    ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
824 
  }
1883 eromero 825 
 
2022 eromero 826 
#if defined(PETSC_USE_COMPLEX)
2682 jroman 827 
  if (d->T) {
828 
    ierr = EPSCleanDenseSchur(d->size_H,0,d->S,d->ldS,d->T,d->ldT,d->eigi-d->cX_in_H,d->pX,d->ldpX,PETSC_TRUE);CHKERRQ(ierr);
829 
  }
2022 eromero 830 
#endif
831 
 
1750 eromero 832 
  PetscFunctionReturn(0);
833 
}
834 
 
835 
 
836 
#undef __FUNCT__
2223 jroman 837 
#define __FUNCT__ "dvd_calcpairs_res_0"
1820 eromero 838 
/* Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
2682 jroman 839 
   the norm associated to the Schur pair, where i = r_s..r_e
1820 eromero 840 
*/
2021 eromero 841 
PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
2605 eromero 842 
                             Vec *R)
1636 slepc 843 
{
2605 eromero 844 
  PetscInt        i;
1636 slepc 845 
  PetscErrorCode  ierr;
2605 eromero 846 
  Vec             *BV = d->BV?d->BV:d->V;
1636 slepc 847 
 
848 
  PetscFunctionBegin;
849 
 
2605 eromero 850 
  for (i=r_s; i<r_e; i++) {
2624 eromero 851 
    /* nX(i) <- ||X(i)|| */
852 
    if (d->correctXnorm) {
853 
      /* R(i) <- V*pX(i) */
854 
      ierr = SlepcUpdateVectorsZ(&R[i-r_s], 0.0, 1.0,
855 
        &d->V[-d->cX_in_H], d->size_V+d->cX_in_H,
856 
        &d->pX[d->ldpX*(i+d->cX_in_H)], d->ldpX, d->size_H, 1); CHKERRQ(ierr);
857 
      ierr = VecNorm(R[i-r_s], NORM_2, &d->nX[i]);CHKERRQ(ierr);
858 
    } else
859 
      d->nX[i] = 1.0;
860 
 
2605 eromero 861 
    /* R(i-r_s) <- AV*pX(i) */
862 
    ierr = SlepcUpdateVectorsZ(&R[i-r_s], 0.0, 1.0,
863 
      &d->AV[-d->cX_in_H], d->size_AV+d->cX_in_H,
864 
      &d->pX[d->ldpX*(i+d->cX_in_H)], d->ldpX, d->size_H, 1); CHKERRQ(ierr);
1636 slepc 865 
 
2605 eromero 866 
    /* R(i-r_s) <- R(i-r_s) - eigr(i)*BV*pX(i) */
867 
    ierr = SlepcUpdateVectorsZ(&R[i-r_s], 1.0, -d->eigr[i+d->cX_in_H],
868 
      &BV[-d->cX_in_H], d->size_V+d->cX_in_H,
869 
      &d->pX[d->ldpX*(i+d->cX_in_H)], d->ldpX, d->size_H, 1); CHKERRQ(ierr);
1636 slepc 870 
  }
871 
 
2021 eromero 872 
  ierr = d->calcpairs_proj_res(d, r_s, r_e, R); CHKERRQ(ierr);
1883 eromero 873 
 
874 
  PetscFunctionReturn(0);
875 
}
876 
 
877 
#undef __FUNCT__
878 
#define __FUNCT__ "dvd_calcpairs_proj_res"
2021 eromero 879 
PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d, PetscInt r_s,
880 
                                      PetscInt r_e, Vec *R)
1883 eromero 881 
{
882 
  PetscInt        i;
883 
  PetscErrorCode  ierr;
2216 jroman 884 
  PetscBool       lindep;
1883 eromero 885 
  Vec             *cX;
886 
 
887 
  PetscFunctionBegin;
888 
 
1829 eromero 889 
  /* If exists the BcX, R <- orth(BcX, R), nR[i] <- ||R[i]|| */
890 
  if (d->BcX)
891 
    cX = d->BcX;
892 
 
1795 eromero 893 
  /* If exists left subspace, R <- orth(cY, R), nR[i] <- ||R[i]|| */
1829 eromero 894 
  else if (d->cY)
895 
    cX = d->cY;
1751 eromero 896 
 
1795 eromero 897 
  /* If fany configurations, R <- orth(cX, R), nR[i] <- ||R[i]|| */
1829 eromero 898 
  else if (!(DVD_IS(d->sEP, DVD_EP_STD) && DVD_IS(d->sEP, DVD_EP_HERMITIAN)))
899 
    cX = d->cX;
900 
 
901 
  /* Otherwise, nR[i] <- ||R[i]|| */
902 
  else
903 
    cX = PETSC_NULL;
904 
 
2779 eromero 905 
  if (cX) {
906 
    Vec auxV;
907 
    ierr = VecDuplicate(d->auxV[0],&auxV);CHKERRQ(ierr);
908 
    for (i=0; i<r_e-r_s; i++) {
909 
      if (cX && d->orthoV_type == DVD_ORTHOV_BOneMV) {
910 
        ierr = IPBOrthogonalize(d->ipV,d->eps->nds,d->eps->DS,d->BDS,d->nBDS,d->size_cX,PETSC_NULL,d->cX,d->real_BV,d->nBcX,R[i],auxV,PETSC_NULL,&d->nR[r_s+i],&lindep);CHKERRQ(ierr);
911 
      } else {
912 
        ierr = IPOrthogonalize(d->ipI, 0, PETSC_NULL, d->size_cX, PETSC_NULL, cX, R[i], PETSC_NULL, &d->nR[r_s+i], &lindep); CHKERRQ(ierr);
913 
      }
914 
      if(lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) {
915 
        ierr = PetscInfo2(d->eps,"The computed eigenvector residual %D is too low, %G!\n",r_s+i,d->nR[r_s+i]);CHKERRQ(ierr);
916 
      }
1795 eromero 917 
    }
2779 eromero 918 
    ierr = VecDestroy(&auxV);CHKERRQ(ierr);
919 
  }
920 
  if (!cX || (cX && d->orthoV_type == DVD_ORTHOV_BOneMV)) {
2682 jroman 921 
    for(i=0; i<r_e-r_s; i++) {
922 
      ierr = VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
923 
    }
924 
    for(i=0; i<r_e-r_s; i++) {
925 
      ierr = VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
926 
    }
1751 eromero 927 
  }
1636 slepc 928 
 
929 
  PetscFunctionReturn(0);
930 
}
931 
 
2682 jroman 932 
#undef __FUNCT__
933 
#define __FUNCT__ "dvd_calcpairs_eig_res_0"
934 
/* Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
935 
   the norm associated to the eigenpair, where i = r_s..r_e
936 
   R, vectors of Vec of size r_e-r_s,
937 
   auxV, PetscMax(r_e+cX_in_H, 2*(r_e-r_s)) vectors,
938 
   auxS, auxiliar vector of size (d->size_cX+r_e)^2+6(d->size_cX+r_e)+(r_e-r_s)*d->size_H
939 
*/
940 
PetscErrorCode dvd_calcpairs_eig_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
941 
{
942 
  PetscInt        i,ldpX,size_in;
943 
  PetscErrorCode  ierr;
944 
  Vec             *Bx;
945 
  PetscScalar     *pX,*pX0;
946 
  DvdReduction    r;
947 
  DvdReductionChunk
948 
                  ops[2];
949 
  DvdMult_copy_func
950 
                  sr[2];
951 
#if !defined(PETSC_USE_COMPLEX)
952 
  PetscScalar     b[8];
953 
  Vec             X[4];
954 
#endif
955 
 
956 
  PetscFunctionBegin;
957 
 
958 
  /* Quick return */
959 
  if (!d->cS) { PetscFunctionReturn(0); }
960 
 
961 
  size_in = (d->size_cX+r_e)*(d->cX_in_AV+r_e)*(d->cT?2:1);
962 
  /* Check consistency */
963 
  if (d->size_auxV < PetscMax(2*(r_e-r_s),d->cX_in_AV+r_e) || d->size_auxS <
964 
      (d->size_cX+r_e)*(d->size_cX+r_e) /* pX */ + PetscMax(PetscMax(
965 
        (d->size_cX+r_e)*6 /* dvd_compute_eigenvectors */,
966 
        d->size_H*(r_e-r_s) /* pX0 */),
967 
        2*size_in /* SlepcAllReduceSum */ )) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
968 
 
969 
  /* Prepare reductions */
970 
  ierr = SlepcAllReduceSumBegin(ops,2,d->auxS,d->auxS+size_in,size_in,&r,((PetscObject)d->V[0])->comm);CHKERRQ(ierr);
971 
 
972 
  /* auxV <- AV * pX(0:r_e+cX_in_H) */
973 
  ierr = SlepcUpdateVectorsZ(d->auxV,0.0,1.0,d->AV-d->cX_in_AV,d->size_AV+d->cX_in_AV,d->pX,d->ldpX,d->size_H,d->cX_in_AV+r_e);CHKERRQ(ierr);
974 
 
975 
  /* cS(:, size_cS:) <- cX' * auxV */
976 
  ierr = VecsMultS(&d->cS[d->ldcS*d->size_cS],0,d->ldcS,d->cY?d->cY:d->cX,0,d->size_cX+r_e,d->auxV,0,d->cX_in_AV+r_e,&r,&sr[0]);CHKERRQ(ierr);
977 
 
978 
  if (d->cT) {
979 
    /* R <- BV * pX(0:r_e+cX_in_H) */
980 
    ierr = SlepcUpdateVectorsZ(d->auxV,0.0,1.0,d->BV-d->cX_in_BV,d->size_BV+d->cX_in_BV,d->pX,d->ldpX,d->size_G,d->cX_in_BV+r_e);CHKERRQ(ierr);
981 
 
982 
    /* cS(:, size_cS:) <- cX' * auxV */
983 
    ierr = VecsMultS(&d->cT[d->ldcT*d->size_cT],0,d->ldcT,d->cY?d->cY:d->cX,0,d->size_cY+r_e,d->auxV,0,d->cX_in_BV+r_e,&r,&sr[1]);CHKERRQ(ierr);
984 
  }
985 
 
986 
  /* Do reductions */
987 
  ierr = SlepcAllReduceSumEnd(&r); CHKERRQ(ierr);
988 
 
989 
  /* qX <- eig(S,T) */
990 
  pX = d->auxS; d->auxS+= (d->size_cX+r_e)*(d->size_cX+r_e); d->size_auxS -= (d->size_cX+r_e)*(d->size_cX+r_e);
991 
  ldpX = d->size_cX+r_e;
992 
  pX0 = d->auxS;
993 
  ierr = EPSCleanDenseSchur(d->size_cX+r_e,0,d->cS,d->ldcS,d->cT,d->ldcT,d->ceigi,pX,ldpX,PETSC_FALSE);CHKERRQ(ierr);
994 
  ierr = dvd_compute_eigenvectors(d->size_cX+r_e,d->cS,d->ldcS,d->cT,d->ldcT,pX,ldpX,PETSC_NULL,0,d->auxS,d->size_auxS,PETSC_TRUE);CHKERRQ(ierr);
995 
 
996 
  /* pX[i] <- pX[i] / ||pX[i]|| */
997 
  pX = &pX[ldpX*d->size_cX+r_s];
998 
  ierr = SlepcDenseNorm(pX,ldpX,d->size_cX+r_e,r_e-r_s,d->eigi+r_s);CHKERRQ(ierr);
999 
 
1000 
  /* pX0 <- d->pX(0:d->cX_in_AV+r_e-1) * pX(size_cX-cX_in_H:) */
1001 
  ierr = SlepcDenseMatProd(pX0,d->size_H,0.0,1.0,d->pX,d->ldpX,d->size_H,d->cX_in_AV+r_e,PETSC_FALSE,&pX[d->size_cX-d->cX_in_H],ldpX,d->cX_in_AV+r_e,r_e-r_s,PETSC_FALSE);CHKERRQ(ierr);
1002 
 
1003 
  /* auxV <- cX(0:size_cX-cX_in_AV)*pX + V*pX0 */
1004 
  ierr = SlepcUpdateVectorsZ(d->auxV,0.0,1.0,d->cX,d->size_cX-d->cX_in_AV,pX,ldpX,d->size_cX-d->cX_in_AV,r_e-r_s);CHKERRQ(ierr);
1005 
  ierr = SlepcUpdateVectorsZ(d->auxV,d->size_cX-d->cX_in_AV==0?0.0:1.0,1.0,d->V-d->cX_in_AV,d->size_V+d->cX_in_AV,pX0,d->size_H,d->size_H,r_e-r_s);CHKERRQ(ierr);
1006 
  d->auxS-= (d->size_cX+r_e)*(d->size_cX+r_e); d->size_auxS += (d->size_cX+r_e)*(d->size_cX+r_e);
1007 
  /* R <- A*auxV */
1008 
  for (i=0; i<r_e-r_s; i++) {
1009 
    ierr = MatMult(d->A,d->auxV[i],R[i]);CHKERRQ(ierr);
1010 
  }
1011 
  /* Bx <- B*auxV */
1012 
  if (d->B) {
1013 
    Bx = &d->auxV[r_e-r_s];
1014 
    for (i=0; i<r_e-r_s; i++) {
1015 
      ierr = MatMult(d->B,d->auxV[i],Bx[i]);CHKERRQ(ierr);
1016 
    }
1017 
  } else {
1018 
    Bx = d->auxV;
1019 
  }
1020 
  /* R <- (A - eig*B)*V*pX */
1021 
  for(i=0; i<r_e-r_s; i++) {
1022 
#if !defined(PETSC_USE_COMPLEX)
1023 
    if (d->eigi[r_s+i] != 0.0) {
1024 
      /* [Ax_i Ax_i+1 Bx_i Bx_i+1]*= [   1        0
1025 
 
1026 
                                      -eigr_i -eigi_i
1027 
                                       eigi_i -eigr_i] */
1028 
      b[0] = b[5] = 1.0;
1029 
      b[2] = b[7] = -d->eigr[r_s+i];
1030 
      b[6] = -(b[3] = d->eigi[r_s+i]);
1031 
      b[1] = b[4] = 0.0;
1032 
      X[0] = R[i]; X[1] = R[i+1]; X[2] = Bx[i]; X[3] = Bx[i+1];
1033 
      ierr = SlepcUpdateVectorsD(X,4,1.0,b,4,4,2,d->auxS,d->size_auxS);CHKERRQ(ierr);
1034 
      i++;
1035 
    } else
1036 
#endif
1037 
    {
1038 
      /* R <- Ax -eig*Bx */
1039 
      ierr = VecAXPBY(R[i], -d->eigr[r_s+i], 1.0, Bx[i]); CHKERRQ(ierr);
1040 
    }
1041 
  }
1042 
 
1043 
  /* nR <- ||R|| */
1044 
  for(i=0; i<r_e-r_s; i++) {
1045 
    ierr = VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]);CHKERRQ(ierr);
1046 
  }
1047 
  for(i=0; i<r_e-r_s; i++) {
1048 
    ierr = VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]);CHKERRQ(ierr);
1049 
  }
1050 
 
1051 
  PetscFunctionReturn(0);
1052 
}
1053 
 
1054 
 
2605 eromero 1055 
/**** Pattern routines ********************************************************/
1056 
 
1057 
/* BV <- BV*MTX */
1735 eromero 1058 
#undef __FUNCT__
2605 eromero 1059 
#define __FUNCT__ "dvd_calcpairs_updateBV0_gen"
1060 
PETSC_STATIC_INLINE PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,Vec *real_BV,PetscInt *size_cBV,Vec **BV,PetscInt *size_BV,PetscInt *max_size_BV,PetscBool BV_shift,PetscInt *cX_in_proj,PetscScalar *MTX)
1734 eromero 1061 
{
1062 
  PetscErrorCode  ierr;
2605 eromero 1063 
  PetscInt        cMT, rm, cp, tra_s, i;
1064 
  Vec             *nBV;
1734 eromero 1065 
 
1066 
  PetscFunctionBegin;
1067 
 
2605 eromero 1068 
  if (!real_BV || !*BV || (d->V_tra_s == 0 && d->V_tra_e == 0)) { PetscFunctionReturn(0); }
1069 
 
1070 
  if (d->V_tra_s > d->max_cX_in_proj && !BV_shift) {
1071 
    tra_s = PetscMax(d->V_tra_s-d->max_cX_in_proj, 0);
1072 
    cMT = d->V_tra_e - tra_s;
1073 
    rm = d->V_tra_s - tra_s;
1074 
    cp = PetscMin(d->max_cX_in_proj - rm, *cX_in_proj);
1075 
    nBV = real_BV+d->max_cX_in_proj;
1076 
    /* BV(-cp-rm:-1-rm) <- BV(-cp:-1) */
1077 
    for (i=-cp; i<0; i++) {
1078 
      ierr = VecCopy((*BV)[i], nBV[i-rm]);CHKERRQ(ierr);
2555 eromero 1079 
    }
2605 eromero 1080 
    /* BV(-rm:) <- BV*MTX(tra_s:V_tra_e-1) */
2624 eromero 1081 
    ierr = SlepcUpdateVectorsZ(&nBV[-rm], 0.0, 1.0, *BV-*cX_in_proj, *size_BV+*cX_in_proj, &MTX[d->ldMTX*tra_s], d->ldMTX, d->size_MT, cMT); CHKERRQ(ierr);
1082 
    *size_BV = d->V_tra_e  - d->V_tra_s;
2605 eromero 1083 
    *max_size_BV-= nBV - *BV;
1084 
    *BV = nBV;
1085 
    if (cX_in_proj && d->max_cX_in_proj>0) *cX_in_proj = cp+rm;
1086 
  } else if (d->V_tra_s <= d->max_cX_in_proj || BV_shift) {
1087 
    /* [BcX BV] <- [BcX BV*MTX] */
1088 
    ierr = SlepcUpdateVectorsZ(*BV-*cX_in_proj, 0.0, 1.0, *BV-*cX_in_proj, *size_BV+*cX_in_proj, MTX, d->ldMTX, d->size_MT, d->V_tra_e); CHKERRQ(ierr);
2624 eromero 1089 
    *BV+= d->V_tra_s-*cX_in_proj;
1090 
    *max_size_BV-= d->V_tra_s-*cX_in_proj;
1091 
    *size_BV = d->V_tra_e  - d->V_tra_s;
1092 
    if (size_cBV && BV_shift) *size_cBV = *BV - real_BV;
1093 
    if (d->max_cX_in_proj>0) *cX_in_proj = PetscMin(*BV - real_BV, d->max_cX_in_proj);
2605 eromero 1094 
  } else { /* !BV_shift */
1095 
    /* BV <- BV*MTX(V_tra_s:) */
1096 
    ierr = SlepcUpdateVectorsZ(*BV, 0.0, 1.0, *BV, *size_BV,
1097 
      &MTX[d->V_tra_s*d->ldMTX], d->ldMTX, d->size_MT, d->V_tra_e-d->V_tra_s);
1795 eromero 1098 
    CHKERRQ(ierr);
2605 eromero 1099 
    *size_BV = d->V_tra_e - d->V_tra_s;
1795 eromero 1100 
  }
1734 eromero 1101 
 
1756 eromero 1102 
  PetscFunctionReturn(0);
1103 
}