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