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