Subversion Repositories slepc-dev

Rev

Go to most recent revision | 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
2116 eromero 15 
   Copyright (c) 2002-2010, Universidad 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);
2021 eromero 38 
PetscErrorCode dvd_calcpairs_updateV(dvdDashboard *d);
39 
PetscErrorCode dvd_calcpairs_updateW(dvdDashboard *d);
40 
PetscErrorCode dvd_calcpairs_updateAV(dvdDashboard *d);
41 
PetscErrorCode dvd_calcpairs_updateBV(dvdDashboard *d);
42 
PetscErrorCode dvd_calcpairs_VtAV_gen(dvdDashboard *d, DvdReduction *r,
43 
                                      DvdMult_copy_func *sr);
44 
PetscErrorCode dvd_calcpairs_VtBV_gen(dvdDashboard *d, DvdReduction *r,
45 
                                      DvdMult_copy_func *sr);
2555 eromero 46 
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d);
2021 eromero 47 
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d);
48 
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d);
49 
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n);
2555 eromero 50 
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n);
2021 eromero 51 
PetscErrorCode dvd_calcpairs_X(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
52 
                               Vec *X);
53 
PetscErrorCode dvd_calcpairs_Y(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
54 
                               Vec *Y);
55 
PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
56 
                                   Vec *R, PetscScalar *auxS, Vec auxV);
57 
PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d, PetscInt r_s,
58 
                                      PetscInt r_e, Vec *R);
59 
PetscErrorCode dvd_calcpairs_updateMatV(Mat A, Vec **AV, PetscInt *size_AV,
2555 eromero 60 
                                        PetscBool doUpdate, PetscBool doNew,
2021 eromero 61 
                                        dvdDashboard *d);
62 
PetscErrorCode dvd_calcpairs_WtMatV_gen(PetscScalar **H, MatType_t sH,
63 
  PetscInt ldH, PetscInt *size_H, PetscScalar *MTY, PetscInt ldMTY,
64 
  PetscScalar *MTX, PetscInt ldMTX, PetscInt rMT, PetscInt cMT, Vec *W,
65 
  Vec *V, PetscInt size_V, PetscScalar *auxS, DvdReduction *r,
66 
  DvdMult_copy_func *sr, dvdDashboard *d);
1619 slepc 67 
 
1735 eromero 68 
/**** Control routines ********************************************************/
1750 eromero 69 
#undef __FUNCT__  
70 
#define __FUNCT__ "dvd_calcpairs_qz"
2021 eromero 71 
PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d, dvdBlackboard *b, IP ipI)
1750 eromero 72 
{
2555 eromero 73 
  PetscErrorCode  ierr;
74 
  PetscInt        i;
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 */
86 
  b->max_nev = PetscMax(b->max_nev, d->nev+1);
87 
#else
88 
  b->max_nev = PetscMax(b->max_nev, d->nev);
89 
#endif
2555 eromero 90 
  b->own_vecs+= b->size_V*(d->B?2:1) + d->eps->nds*(d->ipV_oneMV?1:0);
91 
                                               /* AV, BV?, BDS? */
1757 eromero 92 
  b->own_scalars+= b->max_size_V*b->max_size_V*2*(std_probl?1:2);
93 
                                              /* H, G?, S, T? */
1795 eromero 94 
  b->own_scalars+= b->max_size_V*b->max_size_V*(std_probl?1:2);
1757 eromero 95 
                                              /* pX, pY? */
2555 eromero 96 
  b->own_scalars+= b->max_nev*b->max_nev*(her_probl?0:(std_probl?1:2)); /* cS?, cT?? */
1820 eromero 97 
  b->max_size_auxS = PetscMax(PetscMax(
98 
                              b->max_size_auxS,
99 
                              b->max_size_V*b->max_size_V*4
100 
                                                      /* SlepcReduction */ ),
101 
                              std_probl?0:(b->max_size_V*11+16) /* projeig */);
1750 eromero 102 
 
103 
  /* Setup the step */
104 
  if (b->state >= DVD_STATE_CONF) {
2047 eromero 105 
    d->real_AV = d->AV = b->free_vecs; b->free_vecs+= b->size_V;
106 
    d->max_size_AV = b->size_V;
2555 eromero 107 
    d->max_size_proj = b->max_size_V;
1750 eromero 108 
    d->H = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
109 
    d->real_H = d->H;
110 
    d->pX = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
111 
    d->S = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
2555 eromero 112 
    if (!her_probl) {
113 
      d->cS = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
114 
      d->max_size_cS = b->max_nev;
115 
    } else {
116 
      d->cS = PETSC_NULL;
117 
      d->max_size_cS = 0;
118 
    }
1966 eromero 119 
    d->ldcS = b->max_nev;
1795 eromero 120 
    d->ipV = ipI;
121 
    d->ipW = ipI;
2555 eromero 122 
    if (d->ipV_oneMV) {
123 
      d->BDS = b->free_vecs; b->free_vecs+= d->eps->nds;
124 
      for (i=0; i<d->eps->nds; i++) {
125 
        ierr = MatMult(d->B, d->eps->DS[i], d->BDS[i]); CHKERRQ(ierr);
126 
      }
127 
    }
1883 eromero 128 
    if (d->B) {
2047 eromero 129 
      d->real_BV = d->BV = b->free_vecs; b->free_vecs+= b->size_V;
2555 eromero 130 
    } else {
131 
      d->real_BV = d->BV = PETSC_NULL;
1883 eromero 132 
    }
133 
    if (!std_probl) {
1750 eromero 134 
      d->G = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
135 
      d->real_G = d->G;
136 
      d->T = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
1966 eromero 137 
      d->cT = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
138 
      d->ldcT = b->max_nev;
1795 eromero 139 
      d->pY = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
2247 eromero 140 
    } else {
141 
      d->real_G = d->G = PETSC_NULL;
142 
      d->T = PETSC_NULL;
143 
      d->cT = PETSC_NULL;
144 
      d->ldcT = 0;
145 
      d->pY = PETSC_NULL;
1750 eromero 146 
    }
147 
 
2555 eromero 148 
    d->calcPairs = dvd_calcpairs_proj;
1750 eromero 149 
    d->calcpairs_residual = dvd_calcpairs_res_0;
1883 eromero 150 
    d->calcpairs_proj_res = dvd_calcpairs_proj_res;
2555 eromero 151 
    d->calcpairs_selectPairs = PETSC_NULL;
1750 eromero 152 
    d->calcpairs_X = dvd_calcpairs_X;
1809 eromero 153 
    d->calcpairs_Y = dvd_calcpairs_Y;
1751 eromero 154 
    d->ipI = ipI;
2555 eromero 155 
    d->doNotUpdateBV = PETSC_FALSE;
2244 eromero 156 
    DVD_FL_ADD(d->startList, dvd_calcpairs_qz_start);
1750 eromero 157 
  }
158 
 
159 
  PetscFunctionReturn(0);
160 
}
161 
 
162 
 
2244 eromero 163 
#undef __FUNCT__  
164 
#define __FUNCT__ "dvd_calcpairs_qz_start"
165 
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
166 
{
167 
  PetscBool       std_probl, her_probl;
168 
  PetscInt        i;
169 
 
170 
  PetscFunctionBegin;
171 
 
172 
  std_probl = DVD_IS(d->sEP, DVD_EP_STD)?PETSC_TRUE:PETSC_FALSE;
173 
  her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
174 
 
175 
  d->size_AV = 0;
176 
  d->AV = d->real_AV;
177 
  d->max_size_AV = d->max_size_V;
178 
  d->size_H = 0;
179 
  d->H = d->real_H;
2555 eromero 180 
  d->ldH = d->max_size_proj;
2244 eromero 181 
  for (i=0; i<d->max_size_cS*d->max_size_cS; i++) d->cS[i] = 0.0;
182 
  d->size_cX = d->size_cY = 0;
183 
  d->size_BV = 0;
184 
  if (d->B) {
185 
    d->BV = d->real_BV;
186 
    d->max_size_BV = d->max_size_V;
187 
  } else {
188 
    d->BV = PETSC_NULL;
189 
    d->max_size_BV = 0;
190 
  }
191 
  d->size_G = 0;
2247 eromero 192 
  d->G = d->real_G;
2244 eromero 193 
  if (!std_probl) {
194 
    for (i=0; i<d->max_size_cS*d->max_size_cS; i++) d->cT[i] = 0.0;
195 
    /* If the problem is GHEP without B-orthonormalization, active BcX */
196 
    if(her_probl) d->BcX = d->AV;
197 
 
198 
    /* Else, active the left and right converged invariant subspaces */
199 
    else {d->cY = d->AV; d->BcX = PETSC_NULL; }
200 
  }
201 
 
202 
  PetscFunctionReturn(0);
203 
}
204 
 
205 
 
1750 eromero 206 
#undef __FUNCT__
2555 eromero 207 
#define __FUNCT__ "dvd_calcpairs_proj"
208 
PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
1750 eromero 209 
{
1753 eromero 210 
  PetscErrorCode  ierr;
211 
  DvdReduction    r;
212 
  DvdReductionChunk
1795 eromero 213 
                  ops[2];
1753 eromero 214 
  DvdMult_copy_func
1795 eromero 215 
                  sr[2];
216 
  PetscInt        size_in = 2*d->size_V*d->size_V;
217 
  PetscScalar     *in = d->auxS, *out = in+size_in;
1753 eromero 218 
 
1750 eromero 219 
  PetscFunctionBegin;
220 
 
1753 eromero 221 
  /* Prepare reductions */
1795 eromero 222 
  ierr = SlepcAllReduceSumBegin(ops, 2, in, out, size_in, &r,
1753 eromero 223 
                                ((PetscObject)d->V[0])->comm); CHKERRQ(ierr);
224 
 
2555 eromero 225 
  /* Update AV, BV, W and the projected matrices */
2021 eromero 226 
  ierr = dvd_calcpairs_updateV(d); CHKERRQ(ierr);
227 
  ierr = dvd_calcpairs_updateAV(d); CHKERRQ(ierr);
2555 eromero 228 
  if (!d->W) {
229 
    ierr = dvd_calcpairs_VtAV_gen(d, &r, &sr[0]); CHKERRQ(ierr);
230 
    if (d->BV) { ierr = dvd_calcpairs_updateBV(d); CHKERRQ(ierr); }
231 
  } else {
232 
    if (d->BV) { ierr = dvd_calcpairs_updateBV(d); CHKERRQ(ierr); }
233 
    ierr = dvd_calcpairs_updateW(d); CHKERRQ(ierr);
234 
    ierr = dvd_calcpairs_VtAV_gen(d, &r, &sr[0]); CHKERRQ(ierr);
235 
  }
2021 eromero 236 
  if (DVD_ISNOT(d->sEP, DVD_EP_STD)) {
237 
    ierr = dvd_calcpairs_VtBV_gen(d, &r, &sr[1]); CHKERRQ(ierr);
238 
  }
1753 eromero 239 
 
240 
  /* Do reductions */
241 
  ierr = SlepcAllReduceSumEnd(&r); CHKERRQ(ierr);
242 
 
2555 eromero 243 
  /* Perform the transformation on the projected problem */
244 
  if (d->calcpairs_proj_trans) {
245 
    ierr = d->calcpairs_proj_trans(d); CHKERRQ(ierr);
246 
  }
247 
 
1750 eromero 248 
  if (d->MT_type != DVD_MT_IDENTITY) {
249 
    d->MT_type = DVD_MT_IDENTITY;
250 
//    d->pX_type|= DVD_MAT_IDENTITY;
251 
    d->V_tra_s = d->V_tra_e = 0;
252 
  }
2555 eromero 253 
 
254 
  /* Solve the projected problem */
1750 eromero 255 
  d->pX_type = 0;
2021 eromero 256 
//TODO: uncomment this condition
1750 eromero 257 
//  if(d->V_new_e - d->V_new_s > 0) {
2021 eromero 258 
    if (DVD_IS(d->sEP, DVD_EP_STD)) {
2555 eromero 259 
      if (DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
260 
        ierr = dvd_calcpairs_projeig_eig(d); CHKERRQ(ierr);
261 
      } else {
262 
        ierr = dvd_calcpairs_projeig_qz_std(d); CHKERRQ(ierr);
263 
      }
2021 eromero 264 
    } else {
265 
      ierr = dvd_calcpairs_projeig_qz_gen(d); CHKERRQ(ierr);
266 
    }
1750 eromero 267 
//  }
1735 eromero 268 
  d->V_new_s = d->V_new_e;
1619 slepc 269 
 
1829 eromero 270 
  /* Check consistency */
271 
  if ((d->size_V != d->V_new_e) || (d->size_V != d->size_H) ||
272 
      (d->size_V != d->size_AV) || (DVD_ISNOT(d->sEP, DVD_EP_STD) && (
2555 eromero 273 
      (d->size_V != d->size_G) || (d->BV && d->size_V != d->size_BV) ))) {
2214 jroman 274 
    SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
1829 eromero 275 
  }
276 
 
1619 slepc 277 
  PetscFunctionReturn(0);
278 
}
279 
 
1735 eromero 280 
/**** Basic routines **********************************************************/
281 
 
282 
#undef __FUNCT__
1795 eromero 283 
#define __FUNCT__ "dvd_calcpairs_updateV"
2021 eromero 284 
PetscErrorCode dvd_calcpairs_updateV(dvdDashboard *d)
1795 eromero 285 
{
2021 eromero 286 
  PetscErrorCode  ierr;
1829 eromero 287 
  Vec             *cX = d->BcX? d->BcX : ( (d->cY && !d->W)? d->cY : d->cX );
1795 eromero 288 
 
289 
  PetscFunctionBegin;
290 
 
291 
  /* V <- gs([cX f.V(0:f.V_new_s-1)], f.V(V_new_s:V_new_e-1)) */
2555 eromero 292 
  if (d->ipV_oneMV) {
293 
    ierr = dvd_BorthV(d->ipV, d->eps->DS, d->BDS, d->eps->nds, d->cX, d->real_BV,
294 
                      d->size_cX, d->V, d->BV, d->V_new_s, d->V_new_e,
295 
                      d->auxS, d->auxV[0], d->eps->rand); CHKERRQ(ierr);
296 
  } else {
297 
    ierr = dvd_orthV(d->ipV, d->eps->DS, d->eps->nds, cX, d->size_cX, d->V,
2027 jroman 298 
                   d->V_new_s, d->V_new_e, d->auxS, d->auxV[0], d->eps->rand);
2555 eromero 299 
    CHKERRQ(ierr);
300 
  }
1795 eromero 301 
 
2021 eromero 302 
  PetscFunctionReturn(0);
1795 eromero 303 
}
304 
 
305 
#undef __FUNCT__
306 
#define __FUNCT__ "dvd_calcpairs_updateW"
2021 eromero 307 
PetscErrorCode dvd_calcpairs_updateW(dvdDashboard *d)
1795 eromero 308 
{
2021 eromero 309 
  PetscErrorCode  ierr;
1795 eromero 310 
 
311 
  PetscFunctionBegin;
312 
 
313 
  /* Update W */
2021 eromero 314 
  ierr = d->calcpairs_W(d); CHKERRQ(ierr);
1795 eromero 315 
 
316 
  /* W <- gs([cY f.W(0:f.V_new_s-1)], f.W(V_new_s:V_new_e-1)) */
2021 eromero 317 
  ierr = dvd_orthV(d->ipW, PETSC_NULL, 0, d->cY, d->size_cY, d->W, d->V_new_s,
2027 jroman 318 
                   d->V_new_e, d->auxS, d->auxV[0], d->eps->rand);
2021 eromero 319 
  CHKERRQ(ierr);
1795 eromero 320 
 
2021 eromero 321 
  PetscFunctionReturn(0);
1795 eromero 322 
}
323 
 
324 
 
325 
#undef __FUNCT__
1735 eromero 326 
#define __FUNCT__ "dvd_calcpairs_updateAV"
2021 eromero 327 
PetscErrorCode dvd_calcpairs_updateAV(dvdDashboard *d)
1619 slepc 328 
{
2021 eromero 329 
  PetscErrorCode  ierr;
1619 slepc 330 
 
331 
  PetscFunctionBegin;
332 
 
1735 eromero 333 
  /* f.AV(f.V_tra) = f.AV * f.MT; f.AV(f.V_new) = A*f.V(f.V_new) */
2555 eromero 334 
  ierr = dvd_calcpairs_updateMatV(d->A, &d->AV, &d->size_AV, PETSC_TRUE, PETSC_TRUE, d);
335 
  CHKERRQ(ierr);
1619 slepc 336 
 
2021 eromero 337 
  PetscFunctionReturn(0);
1735 eromero 338 
}
1640 slepc 339 
 
1735 eromero 340 
#undef __FUNCT__
341 
#define __FUNCT__ "dvd_calcpairs_updateBV"
2021 eromero 342 
PetscErrorCode dvd_calcpairs_updateBV(dvdDashboard *d)
1735 eromero 343 
{
2021 eromero 344 
  PetscErrorCode  ierr;
1619 slepc 345 
 
1735 eromero 346 
  PetscFunctionBegin;
1619 slepc 347 
 
1735 eromero 348 
  /* f.BV(f.V_tra) = f.BV * f.MT; f.BV(f.V_new) = B*f.V(f.V_new) */
2555 eromero 349 
  ierr = dvd_calcpairs_updateMatV(d->B, &d->BV, &d->size_BV, !d->doNotUpdateBV, !d->ipV_oneMV, d);
350 
  CHKERRQ(ierr);
351 
  d->doNotUpdateBV = PETSC_FALSE;
1619 slepc 352 
 
2021 eromero 353 
  PetscFunctionReturn(0);
1630 slepc 354 
}
1619 slepc 355 
 
1630 slepc 356 
#undef __FUNCT__
1735 eromero 357 
#define __FUNCT__ "dvd_calcpairs_VtAV_gen"
2021 eromero 358 
PetscErrorCode dvd_calcpairs_VtAV_gen(dvdDashboard *d, DvdReduction *r,
359 
                                      DvdMult_copy_func *sr)
1630 slepc 360 
{
2021 eromero 361 
  PetscInt        ldMTY = d->MTY?d->ldMTY:d->ldMTX;
1753 eromero 362 
  /* WARNING: auxS uses space assigned to r */
1795 eromero 363 
  PetscScalar     *auxS = r->out,
364 
                  *MTY = d->MTY?d->MTY:d->MTX;
365 
  Vec             *W = d->W?d->W:d->V;
2021 eromero 366 
  PetscErrorCode  ierr;
1626 slepc 367 
 
1630 slepc 368 
  PetscFunctionBegin;
1619 slepc 369 
 
1795 eromero 370 
  /* f.H = [f.H(f.V_imm,f.V_imm)        f.V(f.V_imm)'*f.AV(f.V_new);
371 
            f.V(f.V_new)'*f.AV(f.V_imm) f.V(f.V_new)'*f.AV(f.V_new) ] */
1756 eromero 372 
  if (DVD_IS(d->sA,DVD_MAT_HERMITIAN))
1747 eromero 373 
    d->sH = DVD_MAT_HERMITIAN | DVD_MAT_IMPLICIT | DVD_MAT_UTRIANG;
1829 eromero 374 
  if ((d->V_imm_e - d->V_imm_s == 0) && (d->V_tra_e - d->V_tra_s == 0))
375 
    d->size_H = 0;
2021 eromero 376 
  ierr = dvd_calcpairs_WtMatV_gen(&d->H, d->sH, d->ldH, &d->size_H,
1795 eromero 377 
                                  &MTY[ldMTY*d->V_tra_s], ldMTY,
2021 eromero 378 
                                  &d->MTX[d->ldMTX*d->V_tra_s], d->ldMTX,
379 
                                  d->size_MT, d->V_tra_e-d->V_tra_s,
380 
                                  W, d->AV, d->size_V,
381 
                                 auxS, r, sr, d); CHKERRQ(ierr);
1619 slepc 382 
 
2021 eromero 383 
  PetscFunctionReturn(0);
1619 slepc 384 
}
385 
 
1795 eromero 386 
 
1630 slepc 387 
#undef __FUNCT__
1735 eromero 388 
#define __FUNCT__ "dvd_calcpairs_VtBV_gen"
2021 eromero 389 
PetscErrorCode dvd_calcpairs_VtBV_gen(dvdDashboard *d, DvdReduction *r,
390 
                                      DvdMult_copy_func *sr)
1619 slepc 391 
{
2021 eromero 392 
  PetscErrorCode  ierr;
393 
  PetscInt        ldMTY = d->MTY?d->ldMTY:d->ldMTX;
1753 eromero 394 
  /* WARNING: auxS uses space assigned to r */
1795 eromero 395 
  PetscScalar     *auxS = r->out,
396 
                  *MTY = d->MTY?d->MTY:d->MTX;
397 
  Vec             *W = d->W?d->W:d->V;
1619 slepc 398 
 
399 
  PetscFunctionBegin;
400 
 
1795 eromero 401 
  /* f.G = [f.G(f.V_imm,f.V_imm)        f.V(f.V_imm)'*f.BV(f.V_new);
402 
            f.V(f.V_new)'*f.BV(f.V_imm) f.V(f.V_new)'*f.BV(f.V_new) ] */
1756 eromero 403 
  if (DVD_IS(d->sB,DVD_MAT_HERMITIAN))
1747 eromero 404 
    d->sG = DVD_MAT_HERMITIAN | DVD_MAT_IMPLICIT | DVD_MAT_UTRIANG;
1829 eromero 405 
  if ((d->V_imm_e - d->V_imm_s == 0) && (d->V_tra_e - d->V_tra_s == 0))
406 
    d->size_G = 0;
2021 eromero 407 
  ierr = dvd_calcpairs_WtMatV_gen(&d->G, d->sG, d->ldH, &d->size_G,
1795 eromero 408 
                                  &MTY[ldMTY*d->V_tra_s], ldMTY,
2021 eromero 409 
                                  &d->MTX[d->ldMTX*d->V_tra_s], d->ldMTX,
410 
                                  d->size_MT, d->V_tra_e-d->V_tra_s,
411 
                                  W, d->BV?d->BV:d->V, d->size_V,
412 
                                  auxS, r, sr, d); CHKERRQ(ierr);
1630 slepc 413 
 
2021 eromero 414 
  PetscFunctionReturn(0);
1619 slepc 415 
}
416 
 
417 
 
1735 eromero 418 
#undef __FUNCT__
2555 eromero 419 
#define __FUNCT__ "dvd_calcpairs_projeig_eig"
420 
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d)
421 
{
422 
  PetscErrorCode  ierr;
423 
 
424 
  PetscFunctionBegin;
425 
 
426 
  /* S <- H */
427 
  d->ldS = d->ldpX = d->size_H;
428 
  ierr = SlepcDenseCopyTriang(d->S, DVD_MAT_LTRIANG, d->size_H, d->H, d->sH, d->ldH,
429 
                              d->size_H, d->size_H);
430 
 
431 
  /* S = pX' * L * pX */
432 
  ierr = EPSDenseHEP(d->size_H, d->S, d->ldS, d->eigr, d->pX); CHKERRQ(ierr);
433 
 
434 
  d->pX_type = (d->pX_type & !DVD_MAT_IDENTITY) | DVD_MAT_UNITARY;
435 
 
436 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
437 
 
438 
  PetscFunctionReturn(0);
439 
}
440 
 
441 
 
442 
#undef __FUNCT__
1750 eromero 443 
#define __FUNCT__ "dvd_calcpairs_projeig_qz_std"
2021 eromero 444 
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d)
1636 slepc 445 
{
446 
  PetscErrorCode  ierr;
447 
 
448 
  PetscFunctionBegin;
449 
 
1750 eromero 450 
  /* S <- H */
451 
  d->ldS = d->ldpX = d->size_H;
1752 eromero 452 
  ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
453 
                              d->size_H, d->size_H);
1750 eromero 454 
 
455 
  /* S = pX' * H * pX */
2021 eromero 456 
  ierr = EPSDenseHessenberg(d->size_H, 0, d->S, d->ldS, d->pX); CHKERRQ(ierr);
1750 eromero 457 
  ierr = EPSDenseSchur(d->size_H, 0, d->S, d->ldS, d->pX, d->eigr, d->eigi);
458 
  CHKERRQ(ierr);
1636 slepc 459 
 
1750 eromero 460 
  d->pX_type = (d->pX_type & !DVD_MAT_IDENTITY) | DVD_MAT_UNITARY;
1735 eromero 461 
 
2555 eromero 462 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_qz;
463 
 
1636 slepc 464 
  PetscFunctionReturn(0);
465 
}
466 
 
2223 jroman 467 
#undef __FUNCT__
468 
#define __FUNCT__ "dvd_calcpairs_projeig_qz_gen"
1750 eromero 469 
/*
470 
  auxS(dgges) = size_H (beta) + 8*size_H+16 (work)
471 
  auxS(zgges) = size_H (beta) + 1+2*size_H (work) + 8*size_H (rwork)
472 
*/
2021 eromero 473 
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d)
1750 eromero 474 
{
2392 jroman 475 
#if defined(SLEPC_MISSING_LAPACK_GGES)
476 
  PetscFunctionBegin;
477 
  SETERRQ(((PetscObject)(d->eps))->comm,PETSC_ERR_SUP,"GGES - Lapack routine is unavailable.");
478 
#else
1750 eromero 479 
  PetscErrorCode  ierr;
480 
  PetscScalar     *beta = d->auxS;
481 
#if !defined(PETSC_USE_COMPLEX)
482 
  PetscScalar     *auxS = beta + d->size_H;
483 
  PetscBLASInt    n_auxS = d->size_auxS - d->size_H;
484 
#else
1883 eromero 485 
  PetscReal       *auxR = (PetscReal*)(beta + d->size_H);
1750 eromero 486 
  PetscScalar     *auxS = (PetscScalar*)(auxR+8*d->size_H);
1883 eromero 487 
  PetscBLASInt    n_auxS = d->size_auxS - 9*d->size_H;
1750 eromero 488 
#endif
1795 eromero 489 
  PetscInt        i;
1750 eromero 490 
  PetscBLASInt    info,n,a;
491 
 
492 
  PetscFunctionBegin;
493 
  /* S <- H, T <- G */
494 
  d->ldS = d->ldT = d->ldpX = d->ldpY = d->size_H;
1752 eromero 495 
  ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
2139 jroman 496 
                              d->size_H, d->size_H);CHKERRQ(ierr);
1752 eromero 497 
  ierr = SlepcDenseCopyTriang(d->T, 0, d->size_H, d->G, d->sG, d->ldH,
2139 jroman 498 
                              d->size_H, d->size_H);CHKERRQ(ierr);
1750 eromero 499 
 
500 
  /* S = Z'*H*Q, T = Z'*G*Q */
501 
  n = d->size_H;
502 
#if !defined(PETSC_USE_COMPLEX)
1757 eromero 503 
  LAPACKgges_(d->pY?"V":"N", "V", "N", PETSC_NULL, &n, d->S, &n, d->T, &n,
1750 eromero 504 
              &a, d->eigr, d->eigi, beta, d->pY, &n, d->pX, &n,
505 
              auxS, &n_auxS, PETSC_NULL, &info);
506 
#else
1757 eromero 507 
  LAPACKgges_(d->pY?"V":"N", "V", "N", PETSC_NULL, &n, d->S, &n, d->T, &n,
1750 eromero 508 
              &a, d->eigr, beta, d->pY, &n, d->pX, &n,
509 
              auxS, &n_auxS, auxR, PETSC_NULL, &info);
510 
#endif
2214 jroman 511 
  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB, "Error in Lapack GGES %d", info);
1750 eromero 512 
 
513 
  /* eigr[i] <- eigr[i] / beta[i] */
1821 eromero 514 
  for (i=0; i<d->size_H; i++)
515 
    d->eigr[i] /= beta[i],
516 
    d->eigi[i] /= beta[i];
1750 eromero 517 
 
518 
  d->pX_type = (d->pX_type & !DVD_MAT_IDENTITY) | DVD_MAT_UNITARY;
519 
  d->pY_type = (d->pY_type & !DVD_MAT_IDENTITY) | DVD_MAT_UNITARY;
2555 eromero 520 
  d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_qz;
521 
 
1750 eromero 522 
  PetscFunctionReturn(0);
2392 jroman 523 
#endif
1750 eromero 524 
}
525 
 
2555 eromero 526 
#undef __FUNCT__
527 
#define __FUNCT__ "dvd_calcpairs_selectPairs_eig"
528 
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n)
529 
{
530 
  PetscErrorCode  ierr;
1750 eromero 531 
 
2555 eromero 532 
  PetscFunctionBegin;
533 
 
534 
  ierr = EPSSortDenseHEP(d->eps, d->size_H, 0, d->eigr, d->pX, d->ldpX);
535 
  CHKERRQ(ierr);
536 
 
537 
  if (d->calcpairs_eigs_trans) {
538 
    ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
539 
  }
540 
 
541 
  PetscFunctionReturn(0);
542 
}
543 
 
544 
 
1750 eromero 545 
#undef __FUNCT__
546 
#define __FUNCT__ "dvd_calcpairs_selectPairs_qz"
2021 eromero 547 
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n)
1750 eromero 548 
{
1805 eromero 549 
  PetscErrorCode  ierr;
2022 eromero 550 
#if defined(PETSC_USE_COMPLEX)
551 
  PetscScalar     s;
552 
  PetscInt        i, j;
553 
#endif
1750 eromero 554 
  PetscFunctionBegin;
555 
 
1805 eromero 556 
  if ((d->ldpX != d->size_H) ||
557 
      ( d->T &&
558 
        ((d->ldS != d->ldT) || (d->ldpX != d->ldpY) ||
559 
         (d->ldpX != d->size_H)) ) ) {
2214 jroman 560 
     SETERRQ(PETSC_COMM_SELF,1, "Error before ordering eigenpairs");
1750 eromero 561 
  }
562 
 
1829 eromero 563 
  if (d->T) {
2039 eromero 564 
    ierr = EPSSortDenseSchurGeneralized(d->eps, d->size_H, 0, n, d->S, d->T,
1829 eromero 565 
                                        d->ldS, d->pY, d->pX, d->eigr,
566 
                                        d->eigi); CHKERRQ(ierr);
567 
  } else {
568 
    ierr = EPSSortDenseSchur(d->eps, d->size_H, 0, d->S, d->ldS, d->pX,
569 
                             d->eigr, d->eigi); CHKERRQ(ierr);
570 
  }
1750 eromero 571 
 
2021 eromero 572 
  if (d->calcpairs_eigs_trans) {
573 
    ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
574 
  }
1883 eromero 575 
 
2022 eromero 576 
  /* Some functions need the diagonal elements in T be real */
577 
#if defined(PETSC_USE_COMPLEX)
2023 eromero 578 
  if (d->T) for(i=0; i<d->size_H; i++)
579 
    if (PetscImaginaryPart(d->T[d->ldT*i+i]) != 0.0) {
580 
      s = PetscConj(d->T[d->ldT*i+i])/PetscAbsScalar(d->T[d->ldT*i+i]);
581 
      for(j=0; j<=i; j++)
582 
        d->T[d->ldT*i+j] = PetscRealPart(d->T[d->ldT*i+j]*s),
583 
        d->S[d->ldS*i+j]*= s;
584 
      for(j=0; j<d->size_H; j++) d->pX[d->ldpX*i+j]*= s;
585 
    }
2022 eromero 586 
#endif
587 
 
1750 eromero 588 
  PetscFunctionReturn(0);
589 
}
590 
 
591 
 
592 
#undef __FUNCT__
1735 eromero 593 
#define __FUNCT__ "dvd_calcpairs_X"
2021 eromero 594 
PetscErrorCode dvd_calcpairs_X(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
595 
                               Vec *X)
1636 slepc 596 
{
1735 eromero 597 
  PetscInt        i;
598 
  PetscErrorCode  ierr;
1636 slepc 599 
 
600 
  PetscFunctionBegin;
601 
 
1735 eromero 602 
  /* X = V * U(0:n-1) */
1756 eromero 603 
  if (DVD_IS(d->pX_type,DVD_MAT_IDENTITY)) {
1735 eromero 604 
    if (d->V != X) for(i=r_s; i<r_e; i++) {
605 
      ierr = VecCopy(d->V[i], X[i]); CHKERRQ(ierr);
606 
    }
607 
  } else {
1750 eromero 608 
    ierr = SlepcUpdateVectorsZ(X, 0.0, 1.0, d->V, d->size_H, &d->pX[d->ldpX*r_s],
609 
                               d->ldpX, d->size_H, r_e-r_s); CHKERRQ(ierr);
1636 slepc 610 
  }
611 
 
1764 eromero 612 
  /* nX[i] <- ||X[i]|| */
2177 jroman 613 
  if (d->correctXnorm) for(i=0; i<r_e-r_s; i++) {
1764 eromero 614 
    ierr = VecNorm(X[i], NORM_2, &d->nX[r_s+i]); CHKERRQ(ierr);
615 
  } else for(i=0; i<r_e-r_s; i++) {
616 
    d->nX[r_s+i] = 1.0;
617 
  }
618 
 
1636 slepc 619 
  PetscFunctionReturn(0);
620 
}
621 
 
1809 eromero 622 
 
1735 eromero 623 
#undef __FUNCT__
1809 eromero 624 
#define __FUNCT__ "dvd_calcpairs_Y"
2021 eromero 625 
PetscErrorCode dvd_calcpairs_Y(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
626 
                               Vec *Y)
1809 eromero 627 
{
628 
  PetscInt        i, ldpX = d->pY?d->ldpY:d->ldpX;
629 
  PetscErrorCode  ierr;
630 
  Vec             *V = d->W?d->W:d->V;
631 
  PetscScalar     *pX = d->pY?d->pY:d->pX;
632 
 
633 
  PetscFunctionBegin;
634 
 
635 
  /* Y = V * pX(0:n-1) */
636 
  if (DVD_IS(d->pX_type,DVD_MAT_IDENTITY)) {
637 
    if (V != Y) for(i=r_s; i<r_e; i++) {
638 
      ierr = VecCopy(V[i], Y[i]); CHKERRQ(ierr);
639 
    }
640 
  } else {
641 
    ierr = SlepcUpdateVectorsZ(Y, 0.0, 1.0, V, d->size_H, &pX[ldpX*r_s], ldpX,
642 
                               d->size_H, r_e-r_s); CHKERRQ(ierr);
643 
  }
644 
 
645 
  PetscFunctionReturn(0);
646 
}
647 
 
2223 jroman 648 
#undef __FUNCT__
649 
#define __FUNCT__ "dvd_calcpairs_res_0"
1820 eromero 650 
/* Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
651 
   the norm, where
652 
   i <- r_s..r_e,
653 
   UL, auxiliar scalar matrix of size size_H*(r_e-r_s),
654 
   auxV, auxiliar global vector.
655 
*/
2021 eromero 656 
PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
1751 eromero 657 
                             Vec *R, PetscScalar *UL, Vec auxV)
1636 slepc 658 
{
1735 eromero 659 
  PetscInt        i, j;
1636 slepc 660 
  PetscErrorCode  ierr;
661 
 
662 
  PetscFunctionBegin;
663 
 
1764 eromero 664 
  /* If the eigenproblem is not reduced to standard */
1820 eromero 665 
  if ((d->B == PETSC_NULL) || DVD_ISNOT(d->sEP, DVD_EP_STD)) {
1735 eromero 666 
    /* UL = f.U(0:n-1) * diag(f.pL(0:n-1)) */
667 
    for(i=r_s; i<r_e; i++) for(j=0; j<d->size_H; j++)
1820 eromero 668 
      UL[d->size_H*(i-r_s)+j] = d->pX[d->ldpX*i+j]*d->eigr[i];
1636 slepc 669 
 
1883 eromero 670 
    if (d->B == PETSC_NULL) {
1820 eromero 671 
      /* R <- V * UL */
672 
      ierr = SlepcUpdateVectorsZ(R, 0.0, 1.0, d->V, d->size_V, UL, d->size_H,
673 
                                 d->size_H, r_e-r_s); CHKERRQ(ierr);
1636 slepc 674 
    } else {
1820 eromero 675 
      /* R <- BV * UL */
676 
      ierr = SlepcUpdateVectorsZ(R, 0.0, 1.0, d->BV, d->size_BV, UL,
677 
                                 d->size_H, d->size_H, r_e-r_s);
1636 slepc 678 
      CHKERRQ(ierr);
679 
    }
1764 eromero 680 
    /* R <- AV*U - R */
1735 eromero 681 
    ierr = SlepcUpdateVectorsZ(R, -1.0, 1.0, d->AV, d->size_AV,
1750 eromero 682 
                               &d->pX[d->ldpX*r_s], d->ldpX, d->size_H, r_e-r_s);
1735 eromero 683 
    CHKERRQ(ierr);
1764 eromero 684 
 
685 
  /* If the problem was reduced to standard, R[i] = B*X[i] */
686 
  } else {
687 
    /* R[i] <- R[i] * eigr[i] */
688 
    for(i=r_s; i<r_e; i++) {
689 
      ierr = VecScale(R[i-r_s], d->eigr[i]); CHKERRQ(ierr);
690 
    }
691 
 
692 
    /* R <- AV*U - R */
693 
    ierr = SlepcUpdateVectorsZ(R, -1.0, 1.0, d->AV, d->size_AV,
694 
                               &d->pX[d->ldpX*r_s], d->ldpX, d->size_H, r_e-r_s);
695 
    CHKERRQ(ierr);
1636 slepc 696 
  }
697 
 
2021 eromero 698 
  ierr = d->calcpairs_proj_res(d, r_s, r_e, R); CHKERRQ(ierr);
1883 eromero 699 
 
700 
  PetscFunctionReturn(0);
701 
}
702 
 
703 
#undef __FUNCT__
704 
#define __FUNCT__ "dvd_calcpairs_proj_res"
2021 eromero 705 
PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d, PetscInt r_s,
706 
                                      PetscInt r_e, Vec *R)
1883 eromero 707 
{
708 
  PetscInt        i;
709 
  PetscErrorCode  ierr;
2216 jroman 710 
  PetscBool       lindep;
1883 eromero 711 
  Vec             *cX;
712 
 
713 
  PetscFunctionBegin;
714 
 
1829 eromero 715 
  /* If exists the BcX, R <- orth(BcX, R), nR[i] <- ||R[i]|| */
716 
  if (d->BcX)
717 
    cX = d->BcX;
718 
 
1795 eromero 719 
  /* If exists left subspace, R <- orth(cY, R), nR[i] <- ||R[i]|| */
1829 eromero 720 
  else if (d->cY)
721 
    cX = d->cY;
1751 eromero 722 
 
1795 eromero 723 
  /* If fany configurations, R <- orth(cX, R), nR[i] <- ||R[i]|| */
1829 eromero 724 
  else if (!(DVD_IS(d->sEP, DVD_EP_STD) && DVD_IS(d->sEP, DVD_EP_HERMITIAN)))
725 
    cX = d->cX;
726 
 
727 
  /* Otherwise, nR[i] <- ||R[i]|| */
728 
  else
729 
    cX = PETSC_NULL;
730 
 
731 
  if (cX) for (i=0; i<r_e-r_s; i++) {
1805 eromero 732 
    ierr = IPOrthogonalize(d->ipI, 0, PETSC_NULL, d->size_cX, PETSC_NULL,
1829 eromero 733 
                           cX, R[i], PETSC_NULL, &d->nR[r_s+i], &lindep);
1762 eromero 734 
    CHKERRQ(ierr);
2177 jroman 735 
    if(lindep || (d->nR[r_s+i] < PETSC_MACHINE_EPSILON)) {
2499 jroman 736 
      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 737 
    }
1764 eromero 738 
 
1751 eromero 739 
  } else for(i=0; i<r_e-r_s; i++) {
1759 eromero 740 
    ierr = VecNorm(R[i], NORM_2, &d->nR[r_s+i]); CHKERRQ(ierr);
1751 eromero 741 
  }
1636 slepc 742 
 
743 
  PetscFunctionReturn(0);
744 
}
745 
 
1735 eromero 746 
/**** Patterns implementation *************************************************/
747 
#undef __FUNCT__
2223 jroman 748 
#define __FUNCT__ "dvd_calcpairs_updateMatV"
2021 eromero 749 
PetscErrorCode dvd_calcpairs_updateMatV(Mat A, Vec **AV, PetscInt *size_AV,
2555 eromero 750 
                                        PetscBool doUpdate, PetscBool doNew,
2021 eromero 751 
                                        dvdDashboard *d)
1734 eromero 752 
{
1735 eromero 753 
  PetscInt        i;
1734 eromero 754 
  PetscErrorCode  ierr;
755 
 
756 
  PetscFunctionBegin;
757 
 
1735 eromero 758 
  /* f.AV((0:f.V_tra.size)+f.imm.s) = f.AV * f.U(f.V_tra) */
2555 eromero 759 
  if (doUpdate) {
760 
    if (d->MT_type == DVD_MT_pX) {
761 
      ierr = SlepcUpdateVectorsZ(*AV+d->V_imm_e, 0.0, 1.0, *AV, *size_AV,
762 
                                 &d->pX[d->ldpX*d->V_tra_s], d->ldpX,
763 
                                 *size_AV, d->V_tra_e-d->V_tra_s); CHKERRQ(ierr);
764 
    } else if (d->MT_type == DVD_MT_ORTHO) {
765 
      ierr = SlepcUpdateVectorsZ(*AV+d->V_imm_e, 0.0, 1.0, *AV, *size_AV,
766 
                                 &d->MTX[d->ldMTX*d->V_tra_s], d->ldMTX,
767 
                                 *size_AV, d->V_tra_e-d->V_tra_s); CHKERRQ(ierr);
768 
    }
1735 eromero 769 
  }
770 
  *AV = *AV+d->V_imm_s;
1734 eromero 771 
 
1735 eromero 772 
  /* f.AV(f.V_new) = A*f.V(f.V_new) */
1743 eromero 773 
  if (d->V_imm_e-d->V_imm_s + d->V_tra_e-d->V_tra_s != d->V_new_s) {
2214 jroman 774 
    SETERRQ(((PetscObject)A)->comm,1, "Incompatible dimensions");
1743 eromero 775 
  }
776 
 
2555 eromero 777 
  if (doNew) for (i=d->V_new_s; i<d->V_new_e; i++) {
1735 eromero 778 
    ierr = MatMult(A, d->V[i], (*AV)[i]); CHKERRQ(ierr);
1734 eromero 779 
  }
1735 eromero 780 
  *size_AV = d->V_new_e;
1734 eromero 781 
 
782 
  PetscFunctionReturn(0);
783 
}
784 
 
2223 jroman 785 
#undef __FUNCT__
786 
#define __FUNCT__ "dvd_calcpairs_WtMatV_gen"
1735 eromero 787 
/*
1795 eromero 788 
  Compute f.H = [MTY'*H*MTX     W(tra)'*V(new);
789 
                 W(new)'*V(tra) W(new)'*V(new) ]
790 
  where
791 
  tra = 0:cMT-1,
792 
  new = cMT:size_V-1,
1735 eromero 793 
  ldH, the leading dimension of H,
1795 eromero 794 
  auxS, auxiliary scalar vector of size ldH*max(tra,size_V),
1735 eromero 795 
  */
2021 eromero 796 
PetscErrorCode dvd_calcpairs_WtMatV_gen(PetscScalar **H, MatType_t sH,
797 
  PetscInt ldH, PetscInt *size_H, PetscScalar *MTY, PetscInt ldMTY,
798 
  PetscScalar *MTX, PetscInt ldMTX, PetscInt rMT, PetscInt cMT, Vec *W,
799 
  Vec *V, PetscInt size_V, PetscScalar *auxS, DvdReduction *r,
800 
  DvdMult_copy_func *sr, dvdDashboard *d)
1734 eromero 801 
{
802 
  PetscErrorCode  ierr;
803 
 
804 
  PetscFunctionBegin;
805 
 
1795 eromero 806 
  /* H <- MTY^T * (H * MTX) */
807 
  if (cMT > 0) {
1753 eromero 808 
    ierr = SlepcDenseMatProdTriang(auxS, 0, ldH,
1747 eromero 809 
                                   *H, sH, ldH, *size_H, *size_H, PETSC_FALSE,
1795 eromero 810 
                                   MTX, 0, ldMTX, rMT, cMT, PETSC_FALSE);
811 
    CHKERRQ(ierr);
1747 eromero 812 
    ierr = SlepcDenseMatProdTriang(*H, sH, ldH,
1795 eromero 813 
                                   MTY, 0, ldMTY, rMT, cMT, PETSC_TRUE,
1867 eromero 814 
                                   auxS, 0, ldH, *size_H, cMT, PETSC_FALSE);
1747 eromero 815 
    CHKERRQ(ierr);
1795 eromero 816 
    *size_H = cMT;
817 
  }
1734 eromero 818 
 
1795 eromero 819 
  /* H = [H              W(tra)'*W(new);
820 
          W(new)'*V(tra) W(new)'*V(new) ] */
821 
  ierr = VecsMultS(*H, sH, ldH, W, *size_H, size_V, V, *size_H, size_V, r, sr);
1756 eromero 822 
  CHKERRQ(ierr);
1795 eromero 823 
  *size_H = size_V;
1756 eromero 824 
 
825 
  PetscFunctionReturn(0);
826 
}