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
545 dsic.upv.es!jroman 1 
/*
2 
      EPS routines related to the solution process.
1376 slepc 3 
 
4 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672 slepc 5 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2575 eromero 6 
   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
1376 slepc 7 
 
1672 slepc 8 
   This file is part of SLEPc.
9 
 
10 
   SLEPc is free software: you can redistribute it and/or modify it under  the
11 
   terms of version 3 of the GNU Lesser General Public License as published by
12 
   the Free Software Foundation.
13 
 
14 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
15 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
16 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
17 
   more details.
18 
 
19 
   You  should have received a copy of the GNU Lesser General  Public  License
20 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
1376 slepc 21 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
545 dsic.upv.es!jroman 22 
*/
1376 slepc 23 
 
2729 jroman 24 
#include <slepc-private/epsimpl.h>   /*I "slepceps.h" I*/
528 dsic.upv.es!antodo 25 
 
2096 eromero 26 
typedef struct {
2760 jroman 27 
  PetscErrorCode (*which_func)(PetscScalar,PetscScalar,PetscScalar,PetscScalar, PetscInt*,void*);
28 
  void *which_ctx;
29 
  ST st;
2096 eromero 30 
} EPSSortForSTData;
31 
 
2223 jroman 32 
#undef __FUNCT__  
33 
#define __FUNCT__ "EPSSortForSTFunc"
2760 jroman 34 
PetscErrorCode EPSSortForSTFunc(PetscScalar ar,PetscScalar ai,
2331 jroman 35 
                                PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
2096 eromero 36 
{
37 
  EPSSortForSTData *data = (EPSSortForSTData*)ctx;
2317 jroman 38 
  PetscErrorCode   ierr;
2096 eromero 39 
 
40 
  PetscFunctionBegin;
2760 jroman 41 
  ierr = STBackTransform(data->st,1,&ar,&ai);CHKERRQ(ierr);
42 
  ierr = STBackTransform(data->st,1,&br,&bi);CHKERRQ(ierr);
43 
  ierr = (*data->which_func)(ar,ai,br,bi,r,data->which_ctx);CHKERRQ(ierr);
2096 eromero 44 
  PetscFunctionReturn(0);
45 
}
46 
 
528 dsic.upv.es!antodo 47 
#undef __FUNCT__  
48 
#define __FUNCT__ "EPSSolve"
49 
/*@
50 
   EPSSolve - Solves the eigensystem.
51 
 
52 
   Collective on EPS
53 
 
54 
   Input Parameter:
55 
.  eps - eigensolver context obtained from EPSCreate()
56 
 
57 
   Options Database:
58 
+   -eps_view - print information about the solver used
59 
.   -eps_view_binary - save the matrices to the default binary file
60 
-   -eps_plot_eigs - plot computed eigenvalues
61 
 
62 
   Level: beginner
63 
 
64 
.seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
65 
@*/
66 
PetscErrorCode EPSSolve(EPS eps)
67 
{
68 
  PetscErrorCode ierr;
1509 slepc 69 
  PetscInt       i;
528 dsic.upv.es!antodo 70 
  PetscReal      re,im;
2089 jroman 71 
  PetscScalar    dot;
2615 jroman 72 
  PetscBool      flg,isfold,iscayley,viewed=PETSC_FALSE;
528 dsic.upv.es!antodo 73 
  PetscViewer    viewer;
74 
  PetscDraw      draw;
75 
  PetscDrawSP    drawsp;
1654 slepc 76 
  STMatMode      matmode;
1713 antodo 77 
  char           filename[PETSC_MAX_PATH_LEN];
2402 jroman 78 
  char           view[10];
2096 eromero 79 
  EPSSortForSTData data;
2089 jroman 80 
  Mat            A,B;
81 
  KSP            ksp;
82 
  Vec            w,x;
528 dsic.upv.es!antodo 83 
 
84 
  PetscFunctionBegin;
2213 jroman 85 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
528 dsic.upv.es!antodo 86 
 
1713 antodo 87 
  flg = PETSC_FALSE;
2216 jroman 88 
  ierr = PetscOptionsGetBool(((PetscObject)eps)->prefix,"-eps_view_binary",&flg,PETSC_NULL);CHKERRQ(ierr);
528 dsic.upv.es!antodo 89 
  if (flg) {
90 
    ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
1422 slepc 91 
    ierr = MatView(A,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));CHKERRQ(ierr);
92 
    if (B) ierr = MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));CHKERRQ(ierr);
528 dsic.upv.es!antodo 93 
  }
94 
 
2402 jroman 95 
  ierr = PetscOptionsGetString(((PetscObject)eps)->prefix,"-eps_view",view,10,&flg);CHKERRQ(ierr);
96 
  if (flg) {
97 
    ierr = PetscStrcmp(view,"before",&viewed);CHKERRQ(ierr);
98 
    if (viewed){
99 
      PetscViewer viewer;
100 
      ierr = PetscViewerASCIIGetStdout(((PetscObject)eps)->comm,&viewer);CHKERRQ(ierr);
101 
      ierr = EPSView(eps,viewer);CHKERRQ(ierr);
102 
    }
103 
  }
104 
 
528 dsic.upv.es!antodo 105 
  /* reset the convergence flag from the previous solves */
106 
  eps->reason = EPS_CONVERGED_ITERATING;
107 
 
2074 jroman 108 
  /* call setup */
2348 jroman 109 
  if (!eps->setupcalled) { ierr = EPSSetUp(eps);CHKERRQ(ierr); }
528 dsic.upv.es!antodo 110 
  eps->evecsavailable = PETSC_FALSE;
1220 slepc 111 
  eps->nconv = 0;
112 
  eps->its = 0;
113 
  for (i=0;i<eps->ncv;i++) eps->eigr[i]=eps->eigi[i]=eps->errest[i]=0.0;
2313 jroman 114 
  ierr = EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->ncv);CHKERRQ(ierr);
780 dsic.upv.es!jroman 115 
 
528 dsic.upv.es!antodo 116 
  ierr = PetscLogEventBegin(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);
780 dsic.upv.es!jroman 117 
 
2823 jroman 118 
  ierr = PetscObjectTypeCompareAny((PetscObject)eps,&flg,EPSARPACK,EPSBLZPACK,EPSTRLAN,EPSBLOPEX,EPSPRIMME,"");CHKERRQ(ierr);
2161 jroman 119 
  if (!flg) {
2760 jroman 120 
    /* temporarily change eigenvalue comparison function */
121 
    data.which_func = eps->which_func;
122 
    data.which_ctx = eps->which_ctx;
123 
    data.st = eps->OP;
2802 jroman 124 
    eps->which_func = (eps->which==EPS_ALL)? SlepcCompareLargestMagnitude: EPSSortForSTFunc;
125 
    eps->which_ctx = (eps->which==EPS_ALL)? PETSC_NULL: &data;
2161 jroman 126 
  }
2096 eromero 127 
 
1947 jroman 128 
  /* call solver */
129 
  ierr = (*eps->ops->solve)(eps);CHKERRQ(ierr);
780 dsic.upv.es!jroman 130 
 
2161 jroman 131 
  if (!flg) {
2760 jroman 132 
    /* restore comparison function */
133 
    eps->which_func = data.which_func;
134 
    eps->which_ctx = data.which_ctx;
2161 jroman 135 
  }
2096 eromero 136 
 
2809 jroman 137 
  /* finished iteration, truncate the Schur decomposition */
2831 jroman 138 
  if (!(eps->which==EPS_ALL)) {    /* cannot change dimension in spectrum slicing */
2820 carcamgo 139 
    ierr = PSSetDimensions(eps->ps,eps->nconv,0,0);CHKERRQ(ierr);
2821 jroman 140 
    /* change the state to raw so that PSVectors() computes eigenvectors from scratch */
2820 carcamgo 141 
    ierr = PSSetState(eps->ps,PS_STATE_RAW);CHKERRQ(ierr);
142 
  }
2809 jroman 143 
 
1654 slepc 144 
  ierr = STGetMatMode(eps->OP,&matmode);CHKERRQ(ierr);
1940 jroman 145 
  if (matmode == ST_MATMODE_INPLACE && eps->ispositive) {
1930 jroman 146 
    /* Purify eigenvectors before reverting operator */
1654 slepc 147 
    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);    
148 
  }
1029 slepc 149 
  ierr = STPostSolve(eps->OP);CHKERRQ(ierr);
528 dsic.upv.es!antodo 150 
  ierr = PetscLogEventEnd(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);
780 dsic.upv.es!jroman 151 
 
2760 jroman 152 
  if (!eps->reason) SETERRQ(((PetscObject)eps)->comm,1,"Internal error, solver returned without setting converged reason");
528 dsic.upv.es!antodo 153 
 
154 
  /* Map eigenvalues back to the original problem, necessary in some
155 
  * spectral transformations */
2053 jroman 156 
  if (eps->ops->backtransform) {
157 
    ierr = (*eps->ops->backtransform)(eps);CHKERRQ(ierr);
158 
  }
528 dsic.upv.es!antodo 159 
 
889 dsic.upv.es!jroman 160 
  /* Adjust left eigenvectors in generalized problems: y = B^T y */
1947 jroman 161 
  if (eps->isgeneralized && eps->leftvecs) {
889 dsic.upv.es!jroman 162 
    ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);
1422 slepc 163 
    ierr = KSPCreate(((PetscObject)eps)->comm,&ksp);CHKERRQ(ierr);
889 dsic.upv.es!jroman 164 
    ierr = KSPSetOperators(ksp,B,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
165 
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
166 
    ierr = MatGetVecs(B,PETSC_NULL,&w);CHKERRQ(ierr);
167 
    for (i=0;i<eps->nconv;i++) {
168 
      ierr = VecCopy(eps->W[i],w);CHKERRQ(ierr);
169 
      ierr = KSPSolveTranspose(ksp,w,eps->W[i]);CHKERRQ(ierr);
170 
    }
2305 jroman 171 
    ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
172 
    ierr = VecDestroy(&w);CHKERRQ(ierr);
889 dsic.upv.es!jroman 173 
  }
1789 antodo 174 
 
2320 jroman 175 
#if !defined(PETSC_USE_COMPLEX)
1789 antodo 176 
  /* reorder conjugate eigenvalues (positive imaginary first) */
177 
  for (i=0; i<eps->nconv-1; i++) {
178 
    if (eps->eigi[i] != 0) {
179 
      if (eps->eigi[i] < 0) {
180 
        eps->eigi[i] = -eps->eigi[i];
181 
        eps->eigi[i+1] = -eps->eigi[i+1];
182 
        if (!eps->evecsavailable) {
183 
          /* the next correction only works with eigenvectors */
184 
          ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
185 
        }
2330 jroman 186 
        ierr = VecScale(eps->V[i+1],-1.0);CHKERRQ(ierr);
1789 antodo 187 
      }
188 
      i++;
189 
    }
190 
  }
191 
#endif
192 
 
2089 jroman 193 
  /* quick and dirty solution for FOLD: recompute eigenvalues as Rayleigh quotients */
2823 jroman 194 
  ierr = PetscObjectTypeCompare((PetscObject)eps->OP,STFOLD,&isfold);CHKERRQ(ierr);
2089 jroman 195 
  if (isfold) {
196 
    ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
197 
    ierr = MatGetVecs(A,&w,PETSC_NULL);CHKERRQ(ierr);
198 
    if (!eps->evecsavailable) { ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr); }
199 
    for (i=0;i<eps->nconv;i++) {
200 
      x = eps->V[i];
201 
      ierr = MatMult(A,x,w);CHKERRQ(ierr);
202 
      ierr = VecDot(w,x,&eps->eigr[i]);CHKERRQ(ierr);
203 
      if (eps->isgeneralized) {
204 
        ierr = MatMult(B,x,w);CHKERRQ(ierr);
205 
        ierr = VecDot(w,x,&dot);CHKERRQ(ierr);
206 
        eps->eigr[i] /= dot;
207 
      }
208 
    }
2305 jroman 209 
    ierr = VecDestroy(&w);CHKERRQ(ierr);
2089 jroman 210 
  }
211 
 
2615 jroman 212 
  /* In the case of Cayley transform, eigenvectors need to be B-normalized */
2823 jroman 213 
  ierr = PetscObjectTypeCompare((PetscObject)eps->OP,STCAYLEY,&iscayley);CHKERRQ(ierr);
2615 jroman 214 
  if (iscayley && eps->isgeneralized && eps->ishermitian) {
215 
    ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);
2629 carcamgo 216 
    ierr = MatGetVecs(B,PETSC_NULL,&w);CHKERRQ(ierr);
2615 jroman 217 
    if (!eps->evecsavailable) { ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr); }
218 
    for (i=0;i<eps->nconv;i++) {
219 
      x = eps->V[i];
220 
      ierr = MatMult(B,x,w);CHKERRQ(ierr);
221 
      ierr = VecDot(w,x,&dot);CHKERRQ(ierr);
2629 carcamgo 222 
      ierr = VecScale(x,1.0/PetscSqrtScalar(dot));CHKERRQ(ierr);
2615 jroman 223 
    }
224 
    ierr = VecDestroy(&w);CHKERRQ(ierr);
225 
  }
226 
 
528 dsic.upv.es!antodo 227 
  /* sort eigenvalues according to eps->which parameter */
2348 jroman 228 
  ierr = EPSSortEigenvalues(eps,eps->nconv,eps->eigr,eps->eigi,eps->perm);CHKERRQ(ierr);
528 dsic.upv.es!antodo 229 
 
2402 jroman 230 
  if (!viewed) {
231 
    ierr = PetscOptionsGetString(((PetscObject)eps)->prefix,"-eps_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
232 
    if (flg && !PetscPreLoadingOn) {
233 
      ierr = PetscViewerASCIIOpen(((PetscObject)eps)->comm,filename,&viewer);CHKERRQ(ierr);
234 
      ierr = EPSView(eps,viewer);CHKERRQ(ierr);
235 
      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
236 
    }
1713 antodo 237 
  }
528 dsic.upv.es!antodo 238 
 
1713 antodo 239 
  flg = PETSC_FALSE;
2216 jroman 240 
  ierr = PetscOptionsGetBool(((PetscObject)eps)->prefix,"-eps_plot_eigs",&flg,PETSC_NULL);CHKERRQ(ierr);
528 dsic.upv.es!antodo 241 
  if (flg) {
242 
    ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Computed Eigenvalues",
243 
                             PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);CHKERRQ(ierr);
244 
    ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
245 
    ierr = PetscDrawSPCreate(draw,1,&drawsp);CHKERRQ(ierr);
2331 jroman 246 
    for (i=0;i<eps->nconv;i++) {
528 dsic.upv.es!antodo 247 
#if defined(PETSC_USE_COMPLEX)
248 
      re = PetscRealPart(eps->eigr[i]);
249 
      im = PetscImaginaryPart(eps->eigi[i]);
250 
#else
251 
      re = eps->eigr[i];
252 
      im = eps->eigi[i];
253 
#endif
254 
      ierr = PetscDrawSPAddPoint(drawsp,&re,&im);CHKERRQ(ierr);
255 
    }
256 
    ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr);
2305 jroman 257 
    ierr = PetscDrawSPDestroy(&drawsp);CHKERRQ(ierr);
258 
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
528 dsic.upv.es!antodo 259 
  }
260 
 
2080 eromero 261 
  /* Remove the initial subspaces */
262 
  eps->nini = 0;
263 
  eps->ninil = 0;
528 dsic.upv.es!antodo 264 
  PetscFunctionReturn(0);
265 
}
266 
 
267 
#undef __FUNCT__  
268 
#define __FUNCT__ "EPSGetIterationNumber"
269 
/*@
270 
   EPSGetIterationNumber - Gets the current iteration number. If the
271 
   call to EPSSolve() is complete, then it returns the number of iterations
272 
   carried out by the solution method.
273 
 
274 
   Not Collective
275 
 
276 
   Input Parameter:
277 
.  eps - the eigensolver context
278 
 
279 
   Output Parameter:
280 
.  its - number of iterations
281 
 
282 
   Level: intermediate
283 
 
1343 slepc 284 
   Note:
285 
   During the i-th iteration this call returns i-1. If EPSSolve() is
286 
   complete, then parameter "its" contains either the iteration number at
287 
   which convergence was successfully reached, or failure was detected.  
288 
   Call EPSGetConvergedReason() to determine if the solver converged or
289 
   failed and why.
528 dsic.upv.es!antodo 290 
 
1343 slepc 291 
.seealso: EPSGetConvergedReason(), EPSSetTolerances()
528 dsic.upv.es!antodo 292 
@*/
1509 slepc 293 
PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
528 dsic.upv.es!antodo 294 
{
295 
  PetscFunctionBegin;
2213 jroman 296 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
528 dsic.upv.es!antodo 297 
  PetscValidIntPointer(its,2);
298 
  *its = eps->its;
299 
  PetscFunctionReturn(0);
300 
}
301 
 
302 
#undef __FUNCT__  
1209 slepc 303 
#define __FUNCT__ "EPSGetOperationCounters"
528 dsic.upv.es!antodo 304 
/*@
1209 slepc 305 
   EPSGetOperationCounters - Gets the total number of operator applications,
306 
   inner product operations and linear iterations used by the ST object
307 
   during the last EPSSolve() call.
528 dsic.upv.es!antodo 308 
 
309 
   Not Collective
310 
 
311 
   Input Parameter:
312 
.  eps - EPS context
313 
 
314 
   Output Parameter:
1209 slepc 315 
+  ops  - number of operator applications
316 
.  dots - number of inner product operations
317 
-  lits - number of linear iterations
528 dsic.upv.es!antodo 318 
 
319 
   Notes:
320 
   When the eigensolver algorithm invokes STApply() then a linear system
321 
   must be solved (except in the case of standard eigenproblems and shift
322 
   transformation). The number of iterations required in this solve is
323 
   accumulated into a counter whose value is returned by this function.
324 
 
1209 slepc 325 
   These counters are reset to zero at each successive call to EPSSolve().
528 dsic.upv.es!antodo 326 
 
327 
   Level: intermediate
328 
 
329 
@*/
1509 slepc 330 
PetscErrorCode EPSGetOperationCounters(EPS eps,PetscInt* ops,PetscInt* dots,PetscInt* lits)
528 dsic.upv.es!antodo 331 
{
1358 slepc 332 
  PetscErrorCode ierr;
333 
 
528 dsic.upv.es!antodo 334 
  PetscFunctionBegin;
2213 jroman 335 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2371 jroman 336 
  if (!eps->OP) { ierr = EPSGetST(eps,&eps->OP);CHKERRQ(ierr); }
1358 slepc 337 
  ierr = STGetOperationCounters(eps->OP,ops,lits);CHKERRQ(ierr);
338 
  if (dots) {
2381 jroman 339 
    if (!eps->ip) { ierr = EPSGetIP(eps,&eps->ip);CHKERRQ(ierr); }
1358 slepc 340 
    ierr = IPGetOperationCounters(eps->ip,dots);CHKERRQ(ierr);
341 
  }
528 dsic.upv.es!antodo 342 
  PetscFunctionReturn(0);
343 
}
344 
 
345 
#undef __FUNCT__  
346 
#define __FUNCT__ "EPSGetConverged"
347 
/*@
348 
   EPSGetConverged - Gets the number of converged eigenpairs.
349 
 
350 
   Not Collective
351 
 
352 
   Input Parameter:
353 
.  eps - the eigensolver context
354 
 
355 
   Output Parameter:
356 
.  nconv - number of converged eigenpairs
357 
 
358 
   Note:
359 
   This function should be called after EPSSolve() has finished.
360 
 
361 
   Level: beginner
362 
 
1811 jroman 363 
.seealso: EPSSetDimensions(), EPSSolve()
528 dsic.upv.es!antodo 364 
@*/
1509 slepc 365 
PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
528 dsic.upv.es!antodo 366 
{
367 
  PetscFunctionBegin;
2213 jroman 368 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1209 slepc 369 
  PetscValidIntPointer(nconv,2);
370 
  *nconv = eps->nconv;
528 dsic.upv.es!antodo 371 
  PetscFunctionReturn(0);
372 
}
373 
 
374 
 
375 
#undef __FUNCT__  
376 
#define __FUNCT__ "EPSGetConvergedReason"
377 
/*@C
378 
   EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
379 
   stopped.
380 
 
381 
   Not Collective
382 
 
383 
   Input Parameter:
384 
.  eps - the eigensolver context
385 
 
386 
   Output Parameter:
387 
.  reason - negative value indicates diverged, positive value converged
388 
 
389 
   Possible values for reason:
390 
+  EPS_CONVERGED_TOL - converged up to tolerance
391 
.  EPS_DIVERGED_ITS - required more than its to reach convergence
2434 jroman 392 
-  EPS_DIVERGED_BREAKDOWN - generic breakdown in method
528 dsic.upv.es!antodo 393 
 
1811 jroman 394 
   Note:
395 
   Can only be called after the call to EPSSolve() is complete.
396 
 
528 dsic.upv.es!antodo 397 
   Level: intermediate
398 
 
399 
.seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
400 
@*/
401 
PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
402 
{
403 
  PetscFunctionBegin;
2213 jroman 404 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1209 slepc 405 
  PetscValidIntPointer(reason,2);
528 dsic.upv.es!antodo 406 
  *reason = eps->reason;
407 
  PetscFunctionReturn(0);
408 
}
409 
 
410 
#undef __FUNCT__  
411 
#define __FUNCT__ "EPSGetInvariantSubspace"
412 
/*@
761 dsic.upv.es!jroman 413 
   EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
414 
   subspace.
528 dsic.upv.es!antodo 415 
 
1811 jroman 416 
   Not Collective, but vectors are shared by all processors that share the EPS
528 dsic.upv.es!antodo 417 
 
418 
   Input Parameter:
419 
.  eps - the eigensolver context
420 
 
421 
   Output Parameter:
422 
.  v - an array of vectors
423 
 
424 
   Notes:
425 
   This function should be called after EPSSolve() has finished.
426 
 
427 
   The user should provide in v an array of nconv vectors, where nconv is
428 
   the value returned by EPSGetConverged().
429 
 
761 dsic.upv.es!jroman 430 
   The first k vectors returned in v span an invariant subspace associated
431 
   with the first k computed eigenvalues (note that this is not true if the
432 
   k-th eigenvalue is complex and matrix A is real; in this case the first
433 
   k+1 vectors should be used). An invariant subspace X of A satisfies Ax
528 dsic.upv.es!antodo 434 
   in X for all x in X (a similar definition applies for generalized
435 
   eigenproblems).
436 
 
437 
   Level: intermediate
438 
 
1936 jroman 439 
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetInvariantSubspaceLeft()
528 dsic.upv.es!antodo 440 
@*/
2331 jroman 441 
PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec *v)
528 dsic.upv.es!antodo 442 
{
443 
  PetscErrorCode ierr;
1509 slepc 444 
  PetscInt       i;
528 dsic.upv.es!antodo 445 
 
446 
  PetscFunctionBegin;
2213 jroman 447 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
796 dsic.upv.es!antodo 448 
  PetscValidPointer(v,2);
2213 jroman 449 
  PetscValidHeaderSpecific(*v,VEC_CLASSID,2);
2762 jroman 450 
  if (!eps->V) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
451 
  if (!eps->ishermitian && eps->evecsavailable) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector,EPSComputeRelativeError or EPSComputeResidualNorm");
1940 jroman 452 
  if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
1804 jroman 453 
    for (i=0;i<eps->nconv;i++) {
454 
      ierr = VecPointwiseDivide(v[i],eps->V[i],eps->D);CHKERRQ(ierr);
455 
      ierr = VecNormalize(v[i],PETSC_NULL);CHKERRQ(ierr);
456 
    }
528 dsic.upv.es!antodo 457 
  }
1804 jroman 458 
  else {
459 
    for (i=0;i<eps->nconv;i++) {
460 
      ierr = VecCopy(eps->V[i],v[i]);CHKERRQ(ierr);
461 
    }
462 
  }
528 dsic.upv.es!antodo 463 
  PetscFunctionReturn(0);
464 
}
465 
 
466 
#undef __FUNCT__  
1936 jroman 467 
#define __FUNCT__ "EPSGetInvariantSubspaceLeft"
780 dsic.upv.es!jroman 468 
/*@
1936 jroman 469 
   EPSGetInvariantSubspaceLeft - Gets an orthonormal basis of the computed left
780 dsic.upv.es!jroman 470 
   invariant subspace (only available in two-sided eigensolvers).
471 
 
1811 jroman 472 
   Not Collective, but vectors are shared by all processors that share the EPS
780 dsic.upv.es!jroman 473 
 
474 
   Input Parameter:
475 
.  eps - the eigensolver context
476 
 
477 
   Output Parameter:
478 
.  v - an array of vectors
479 
 
480 
   Notes:
481 
   This function should be called after EPSSolve() has finished.
482 
 
483 
   The user should provide in v an array of nconv vectors, where nconv is
484 
   the value returned by EPSGetConverged().
485 
 
486 
   The first k vectors returned in v span a left invariant subspace associated
487 
   with the first k computed eigenvalues (note that this is not true if the
488 
   k-th eigenvalue is complex and matrix A is real; in this case the first
489 
   k+1 vectors should be used). A left invariant subspace Y of A satisfies y'A
490 
   in Y for all y in Y (a similar definition applies for generalized
491 
   eigenproblems).
492 
 
493 
   Level: intermediate
494 
 
495 
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetInvariantSubspace
496 
@*/
2326 jroman 497 
PetscErrorCode EPSGetInvariantSubspaceLeft(EPS eps,Vec *v)
780 dsic.upv.es!jroman 498 
{
499 
  PetscErrorCode ierr;
1509 slepc 500 
  PetscInt       i;
780 dsic.upv.es!jroman 501 
 
502 
  PetscFunctionBegin;
2213 jroman 503 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
819 dsic.upv.es!jroman 504 
  PetscValidPointer(v,2);
2213 jroman 505 
  PetscValidHeaderSpecific(*v,VEC_CLASSID,2);
2762 jroman 506 
  if (!eps->leftvecs) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetLeftVectorsWanted");
507 
  if (!eps->W) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
508 
  if (!eps->ishermitian && eps->evecsavailable) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspaceLeft must be called before EPSGetEigenpairLeft,EPSComputeRelativeErrorLeft or EPSComputeResidualNormLeft");
780 dsic.upv.es!jroman 509 
  for (i=0;i<eps->nconv;i++) {
510 
    ierr = VecCopy(eps->W[i],v[i]);CHKERRQ(ierr);
511 
  }
512 
  PetscFunctionReturn(0);
513 
}
514 
 
515 
#undef __FUNCT__  
528 dsic.upv.es!antodo 516 
#define __FUNCT__ "EPSGetEigenpair"
517 
/*@
780 dsic.upv.es!jroman 518 
   EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
519 
   EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
528 dsic.upv.es!antodo 520 
 
1811 jroman 521 
   Not Collective, but vectors are shared by all processors that share the EPS
528 dsic.upv.es!antodo 522 
 
523 
   Input Parameters:
524 
+  eps - eigensolver context
525 
-  i   - index of the solution
526 
 
527 
   Output Parameters:
528 
+  eigr - real part of eigenvalue
529 
.  eigi - imaginary part of eigenvalue
530 
.  Vr   - real part of eigenvector
531 
-  Vi   - imaginary part of eigenvector
532 
 
533 
   Notes:
1389 slepc 534 
   If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
535 
   configured with complex scalars the eigenvalue is stored
761 dsic.upv.es!jroman 536 
   directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
528 dsic.upv.es!antodo 537 
   set to zero).
538 
 
1267 slepc 539 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
528 dsic.upv.es!antodo 540 
   Eigenpairs are indexed according to the ordering criterion established
541 
   with EPSSetWhichEigenpairs().
542 
 
1772 antodo 543 
   The 2-norm of the eigenvector is one unless the problem is generalized
544 
   Hermitian. In this case the eigenvector is normalized with respect to the
545 
   norm defined by the B matrix.
546 
 
528 dsic.upv.es!antodo 547 
   Level: beginner
548 
 
1936 jroman 549 
.seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetEigenvectorLeft(), EPSSolve(),
780 dsic.upv.es!jroman 550 
          EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
528 dsic.upv.es!antodo 551 
@*/
2326 jroman 552 
PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
528 dsic.upv.es!antodo 553 
{
554 
  PetscErrorCode ierr;
780 dsic.upv.es!jroman 555 
 
556 
  PetscFunctionBegin;
2213 jroman 557 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2762 jroman 558 
  if (!eps->eigr || !eps->eigi || !eps->V) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
559 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1936 jroman 560 
  ierr = EPSGetEigenvalue(eps,i,eigr,eigi);CHKERRQ(ierr);
561 
  ierr = EPSGetEigenvector(eps,i,Vr,Vi);CHKERRQ(ierr);
780 dsic.upv.es!jroman 562 
  PetscFunctionReturn(0);
563 
}
564 
 
565 
#undef __FUNCT__  
1936 jroman 566 
#define __FUNCT__ "EPSGetEigenvalue"
780 dsic.upv.es!jroman 567 
/*@
1936 jroman 568 
   EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
780 dsic.upv.es!jroman 569 
 
570 
   Not Collective
571 
 
572 
   Input Parameters:
573 
+  eps - eigensolver context
574 
-  i   - index of the solution
575 
 
576 
   Output Parameters:
577 
+  eigr - real part of eigenvalue
578 
-  eigi - imaginary part of eigenvalue
579 
 
580 
   Notes:
1389 slepc 581 
   If the eigenvalue is real, then eigi is set to zero. If PETSc is
582 
   configured with complex scalars the eigenvalue is stored
780 dsic.upv.es!jroman 583 
   directly in eigr (eigi is set to zero).
584 
 
1267 slepc 585 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
780 dsic.upv.es!jroman 586 
   Eigenpairs are indexed according to the ordering criterion established
587 
   with EPSSetWhichEigenpairs().
588 
 
589 
   Level: beginner
590 
 
591 
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
592 
          EPSGetEigenpair()
593 
@*/
2326 jroman 594 
PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
780 dsic.upv.es!jroman 595 
{
1509 slepc 596 
  PetscInt       k;
780 dsic.upv.es!jroman 597 
 
598 
  PetscFunctionBegin;
2213 jroman 599 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2762 jroman 600 
  if (!eps->eigr || !eps->eigi) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
601 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
780 dsic.upv.es!jroman 602 
  if (!eps->perm) k = i;
603 
  else k = eps->perm[i];
2320 jroman 604 
#if defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 605 
  if (eigr) *eigr = eps->eigr[k];
606 
  if (eigi) *eigi = 0;
607 
#else
608 
  if (eigr) *eigr = eps->eigr[k];
609 
  if (eigi) *eigi = eps->eigi[k];
610 
#endif
611 
  PetscFunctionReturn(0);
612 
}
613 
 
614 
#undef __FUNCT__  
1936 jroman 615 
#define __FUNCT__ "EPSGetEigenvector"
780 dsic.upv.es!jroman 616 
/*@
1936 jroman 617 
   EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
780 dsic.upv.es!jroman 618 
 
1811 jroman 619 
   Not Collective, but vectors are shared by all processors that share the EPS
780 dsic.upv.es!jroman 620 
 
621 
   Input Parameters:
622 
+  eps - eigensolver context
623 
-  i   - index of the solution
624 
 
625 
   Output Parameters:
626 
+  Vr   - real part of eigenvector
627 
-  Vi   - imaginary part of eigenvector
628 
 
629 
   Notes:
1389 slepc 630 
   If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
631 
   configured with complex scalars the eigenvector is stored
780 dsic.upv.es!jroman 632 
   directly in Vr (Vi is set to zero).
633 
 
1267 slepc 634 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
780 dsic.upv.es!jroman 635 
   Eigenpairs are indexed according to the ordering criterion established
636 
   with EPSSetWhichEigenpairs().
637 
 
1772 antodo 638 
   The 2-norm of the eigenvector is one unless the problem is generalized
639 
   Hermitian. In this case the eigenvector is normalized with respect to the
640 
   norm defined by the B matrix.
641 
 
780 dsic.upv.es!jroman 642 
   Level: beginner
643 
 
644 
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
1936 jroman 645 
          EPSGetEigenpair(), EPSGetEigenvectorLeft()
780 dsic.upv.es!jroman 646 
@*/
2326 jroman 647 
PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
780 dsic.upv.es!jroman 648 
{
649 
  PetscErrorCode ierr;
1509 slepc 650 
  PetscInt       k;
528 dsic.upv.es!antodo 651 
 
652 
  PetscFunctionBegin;
2213 jroman 653 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 654 
  if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Vr,3); }
655 
  if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Vi,4); }
2762 jroman 656 
  if (!eps->V) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
657 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
2331 jroman 658 
  if (!eps->evecsavailable && (Vr || Vi)) {
528 dsic.upv.es!antodo 659 
    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
660 
  }  
1789 antodo 661 
  if (!eps->perm) k = i;
662 
  else k = eps->perm[i];
2320 jroman 663 
#if defined(PETSC_USE_COMPLEX)
2331 jroman 664 
  if (Vr) { ierr = VecCopy(eps->V[k],Vr);CHKERRQ(ierr); }
2330 jroman 665 
  if (Vi) { ierr = VecSet(Vi,0.0);CHKERRQ(ierr); }
1789 antodo 666 
#else
667 
  if (eps->eigi[k] > 0) { /* first value of conjugate pair */
2331 jroman 668 
    if (Vr) { ierr = VecCopy(eps->V[k],Vr);CHKERRQ(ierr); }
669 
    if (Vi) { ierr = VecCopy(eps->V[k+1],Vi);CHKERRQ(ierr); }
1789 antodo 670 
  } else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
2331 jroman 671 
    if (Vr) { ierr = VecCopy(eps->V[k-1],Vr);CHKERRQ(ierr); }
1789 antodo 672 
    if (Vi) {
2331 jroman 673 
      ierr = VecCopy(eps->V[k],Vi);CHKERRQ(ierr);
2330 jroman 674 
      ierr = VecScale(Vi,-1.0);CHKERRQ(ierr);
1789 antodo 675 
    }
676 
  } else { /* real eigenvalue */
2331 jroman 677 
    if (Vr) { ierr = VecCopy(eps->V[k],Vr);CHKERRQ(ierr); }
2330 jroman 678 
    if (Vi) { ierr = VecSet(Vi,0.0);CHKERRQ(ierr); }
528 dsic.upv.es!antodo 679 
  }
680 
#endif
681 
  PetscFunctionReturn(0);
682 
}
683 
 
684 
#undef __FUNCT__  
1936 jroman 685 
#define __FUNCT__ "EPSGetEigenvectorLeft"
780 dsic.upv.es!jroman 686 
/*@
1936 jroman 687 
   EPSGetEigenvectorLeft - Gets the i-th left eigenvector as computed by EPSSolve()
780 dsic.upv.es!jroman 688 
   (only available in two-sided eigensolvers).
689 
 
1811 jroman 690 
   Not Collective, but vectors are shared by all processors that share the EPS
780 dsic.upv.es!jroman 691 
 
692 
   Input Parameters:
693 
+  eps - eigensolver context
694 
-  i   - index of the solution
695 
 
696 
   Output Parameters:
697 
+  Wr   - real part of eigenvector
698 
-  Wi   - imaginary part of eigenvector
699 
 
700 
   Notes:
1389 slepc 701 
   If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
702 
   configured with complex scalars the eigenvector is stored
780 dsic.upv.es!jroman 703 
   directly in Wr (Wi is set to zero).
704 
 
1267 slepc 705 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
780 dsic.upv.es!jroman 706 
   Eigenpairs are indexed according to the ordering criterion established
707 
   with EPSSetWhichEigenpairs().
708 
 
709 
   Level: beginner
710 
 
711 
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
1936 jroman 712 
          EPSGetEigenpair(), EPSGetEigenvector()
780 dsic.upv.es!jroman 713 
@*/
2326 jroman 714 
PetscErrorCode EPSGetEigenvectorLeft(EPS eps,PetscInt i,Vec Wr,Vec Wi)
780 dsic.upv.es!jroman 715 
{
716 
  PetscErrorCode ierr;
1509 slepc 717 
  PetscInt       k;
780 dsic.upv.es!jroman 718 
 
719 
  PetscFunctionBegin;
2213 jroman 720 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 721 
  if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Wr,3); }
722 
  if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Wi,4); }
2762 jroman 723 
  if (!eps->leftvecs) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetLeftVectorsWanted");
724 
  if (!eps->W) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
725 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
2331 jroman 726 
  if (!eps->evecsavailable && (Wr || Wi)) {
780 dsic.upv.es!jroman 727 
    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
728 
  }  
1789 antodo 729 
  if (!eps->perm) k = i;
730 
  else k = eps->perm[i];
2320 jroman 731 
#if defined(PETSC_USE_COMPLEX)
2331 jroman 732 
  if (Wr) { ierr = VecCopy(eps->W[k],Wr);CHKERRQ(ierr); }
2330 jroman 733 
  if (Wi) { ierr = VecSet(Wi,0.0);CHKERRQ(ierr); }
1789 antodo 734 
#else
735 
  if (eps->eigi[k] > 0) { /* first value of conjugate pair */
2331 jroman 736 
    if (Wr) { ierr = VecCopy(eps->W[k],Wr);CHKERRQ(ierr); }
737 
    if (Wi) { ierr = VecCopy(eps->W[k+1],Wi);CHKERRQ(ierr); }
1789 antodo 738 
  } else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
2331 jroman 739 
    if (Wr) { ierr = VecCopy(eps->W[k-1],Wr);CHKERRQ(ierr); }
1789 antodo 740 
    if (Wi) {
2331 jroman 741 
      ierr = VecCopy(eps->W[k],Wi);CHKERRQ(ierr);
2330 jroman 742 
      ierr = VecScale(Wi,-1.0);CHKERRQ(ierr);
1789 antodo 743 
    }
744 
  } else { /* real eigenvalue */
2331 jroman 745 
    if (Wr) { ierr = VecCopy(eps->W[k],Wr);CHKERRQ(ierr); }
2330 jroman 746 
    if (Wi) { ierr = VecSet(Wi,0.0);CHKERRQ(ierr); }
780 dsic.upv.es!jroman 747 
  }
748 
#endif
749 
  PetscFunctionReturn(0);
750 
}
751 
 
752 
#undef __FUNCT__  
528 dsic.upv.es!antodo 753 
#define __FUNCT__ "EPSGetErrorEstimate"
754 
/*@
761 dsic.upv.es!jroman 755 
   EPSGetErrorEstimate - Returns the error estimate associated to the i-th
756 
   computed eigenpair.
528 dsic.upv.es!antodo 757 
 
758 
   Not Collective
759 
 
760 
   Input Parameter:
761 
+  eps - eigensolver context
762 
-  i   - index of eigenpair
763 
 
764 
   Output Parameter:
765 
.  errest - the error estimate
766 
 
761 dsic.upv.es!jroman 767 
   Notes:
768 
   This is the error estimate used internally by the eigensolver. The actual
1811 jroman 769 
   error bound can be computed with EPSComputeRelativeError(). See also the users
761 dsic.upv.es!jroman 770 
   manual for details.
771 
 
528 dsic.upv.es!antodo 772 
   Level: advanced
773 
 
774 
.seealso: EPSComputeRelativeError()
775 
@*/
2326 jroman 776 
PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
528 dsic.upv.es!antodo 777 
{
778 
  PetscFunctionBegin;
2213 jroman 779 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 780 
  PetscValidPointer(errest,3);
2762 jroman 781 
  if (!eps->eigr || !eps->eigi) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
782 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
528 dsic.upv.es!antodo 783 
  if (eps->perm) i = eps->perm[i];  
784 
  if (errest) *errest = eps->errest[i];
785 
  PetscFunctionReturn(0);
786 
}
787 
 
780 dsic.upv.es!jroman 788 
#undef __FUNCT__  
789 
#define __FUNCT__ "EPSGetErrorEstimateLeft"
790 
/*@
791 
   EPSGetErrorEstimateLeft - Returns the left error estimate associated to the i-th
792 
   computed eigenpair (only available in two-sided eigensolvers).
528 dsic.upv.es!antodo 793 
 
780 dsic.upv.es!jroman 794 
   Not Collective
795 
 
796 
   Input Parameter:
797 
+  eps - eigensolver context
798 
-  i   - index of eigenpair
799 
 
800 
   Output Parameter:
801 
.  errest - the left error estimate
802 
 
803 
   Notes:
804 
   This is the error estimate used internally by the eigensolver. The actual
1811 jroman 805 
   error bound can be computed with EPSComputeRelativeErrorLeft(). See also the users
780 dsic.upv.es!jroman 806 
   manual for details.
807 
 
808 
   Level: advanced
809 
 
810 
.seealso: EPSComputeRelativeErrorLeft()
811 
@*/
2331 jroman 812 
PetscErrorCode EPSGetErrorEstimateLeft(EPS eps,PetscInt i,PetscReal *errest)
780 dsic.upv.es!jroman 813 
{
814 
  PetscFunctionBegin;
2213 jroman 815 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 816 
  PetscValidPointer(errest,3);
2762 jroman 817 
  if (!eps->eigr || !eps->eigi) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"EPSSolve must be called first");
818 
  if (!eps->leftvecs) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetLeftVectorsWanted");
819 
  if (i<0 || i>=eps->nconv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
780 dsic.upv.es!jroman 820 
  if (eps->perm) i = eps->perm[i];  
821 
  if (errest) *errest = eps->errest_left[i];
822 
  PetscFunctionReturn(0);
823 
}
824 
 
528 dsic.upv.es!antodo 825 
#undef __FUNCT__  
1812 antodo 826 
#define __FUNCT__ "EPSComputeResidualNorm_Private"
827 
/*
828 
   EPSComputeResidualNorm_Private - Computes the norm of the residual vector
829 
   associated with an eigenpair.
830 
*/
2331 jroman 831 
PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,PetscReal *norm)
528 dsic.upv.es!antodo 832 
{
833 
  PetscErrorCode ierr;
2331 jroman 834 
  Vec            u,w;
835 
  Mat            A,B;
2320 jroman 836 
#if !defined(PETSC_USE_COMPLEX)
1812 antodo 837 
  Vec            v;
2331 jroman 838 
  PetscReal      ni,nr;
528 dsic.upv.es!antodo 839 
#endif
840 
 
841 
  PetscFunctionBegin;
842 
  ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
1933 jroman 843 
  ierr = VecDuplicate(eps->V[0],&u);CHKERRQ(ierr);
844 
  ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
1812 antodo 845 
 
2320 jroman 846 
#if !defined(PETSC_USE_COMPLEX)
528 dsic.upv.es!antodo 847 
  if (ki == 0 ||
848 
    PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
849 
#endif
1893 jroman 850 
    ierr = MatMult(A,xr,u);CHKERRQ(ierr);                             /* u=A*x */
528 dsic.upv.es!antodo 851 
    if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
1893 jroman 852 
      if (eps->isgeneralized) { ierr = MatMult(B,xr,w);CHKERRQ(ierr); }
853 
      else { ierr = VecCopy(xr,w);CHKERRQ(ierr); }                    /* w=B*x */
854 
      ierr = VecAXPY(u,-kr,w);CHKERRQ(ierr);                          /* u=A*x-k*B*x */
528 dsic.upv.es!antodo 855 
    }
1893 jroman 856 
    ierr = VecNorm(u,NORM_2,norm);CHKERRQ(ierr);  
2320 jroman 857 
#if !defined(PETSC_USE_COMPLEX)
528 dsic.upv.es!antodo 858 
  } else {
2330 jroman 859 
    ierr = VecDuplicate(eps->V[0],&v);CHKERRQ(ierr);
1893 jroman 860 
    ierr = MatMult(A,xr,u);CHKERRQ(ierr);                             /* u=A*xr */
861 
    if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
862 
      if (eps->isgeneralized) { ierr = MatMult(B,xr,v);CHKERRQ(ierr); }
863 
      else { ierr = VecCopy(xr,v);CHKERRQ(ierr); }                    /* v=B*xr */
864 
      ierr = VecAXPY(u,-kr,v);CHKERRQ(ierr);                          /* u=A*xr-kr*B*xr */
865 
      if (eps->isgeneralized) { ierr = MatMult(B,xi,w);CHKERRQ(ierr); }
866 
      else { ierr = VecCopy(xi,w);CHKERRQ(ierr); }                    /* w=B*xi */
867 
      ierr = VecAXPY(u,ki,w);CHKERRQ(ierr);                           /* u=A*xr-kr*B*xr+ki*B*xi */
868 
    }
869 
    ierr = VecNorm(u,NORM_2,&nr);CHKERRQ(ierr);
870 
    ierr = MatMult(A,xi,u);CHKERRQ(ierr);                             /* u=A*xi */
871 
    if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
872 
      ierr = VecAXPY(u,-kr,w);CHKERRQ(ierr);                          /* u=A*xi-kr*B*xi */
873 
      ierr = VecAXPY(u,-ki,v);CHKERRQ(ierr);                          /* u=A*xi-kr*B*xi-ki*B*xr */
874 
    }
875 
    ierr = VecNorm(u,NORM_2,&ni);CHKERRQ(ierr);
876 
    *norm = SlepcAbsEigenvalue(nr,ni);
2305 jroman 877 
    ierr = VecDestroy(&v);CHKERRQ(ierr);
528 dsic.upv.es!antodo 878 
  }
879 
#endif
880 
 
2305 jroman 881 
  ierr = VecDestroy(&w);CHKERRQ(ierr);
882 
  ierr = VecDestroy(&u);CHKERRQ(ierr);
1812 antodo 883 
  PetscFunctionReturn(0);
884 
}
885 
 
886 
#undef __FUNCT__  
887 
#define __FUNCT__ "EPSComputeResidualNorm"
888 
/*@
889 
   EPSComputeResidualNorm - Computes the norm of the residual vector associated with
890 
   the i-th computed eigenpair.
891 
 
892 
   Collective on EPS
893 
 
894 
   Input Parameter:
895 
.  eps - the eigensolver context
896 
.  i   - the solution index
897 
 
898 
   Output Parameter:
899 
.  norm - the residual norm, computed as ||Ax-kBx||_2 where k is the
900 
   eigenvalue and x is the eigenvector.
901 
   If k=0 then the residual norm is computed as ||Ax||_2.
902 
 
903 
   Notes:
904 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
905 
   Eigenpairs are indexed according to the ordering criterion established
906 
   with EPSSetWhichEigenpairs().
907 
 
908 
   Level: beginner
909 
 
910 
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
911 
@*/
2331 jroman 912 
PetscErrorCode EPSComputeResidualNorm(EPS eps,PetscInt i,PetscReal *norm)
1812 antodo 913 
{
914 
  PetscErrorCode ierr;
2331 jroman 915 
  Vec            xr,xi;
916 
  PetscScalar    kr,ki;
1812 antodo 917 
 
918 
  PetscFunctionBegin;
2213 jroman 919 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2328 jroman 920 
  PetscValidLogicalCollectiveInt(eps,i,2);
1812 antodo 921 
  PetscValidPointer(norm,3);
2330 jroman 922 
  ierr = VecDuplicate(eps->V[0],&xr);CHKERRQ(ierr);
923 
  ierr = VecDuplicate(eps->V[0],&xi);CHKERRQ(ierr);
924 
  ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
925 
  ierr = EPSComputeResidualNorm_Private(eps,kr,ki,xr,xi,norm);CHKERRQ(ierr);
926 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
927 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
528 dsic.upv.es!antodo 928 
  PetscFunctionReturn(0);
929 
}
930 
 
931 
#undef __FUNCT__  
780 dsic.upv.es!jroman 932 
#define __FUNCT__ "EPSComputeResidualNormLeft"
933 
/*@
794 dsic.upv.es!antodo 934 
   EPSComputeResidualNormLeft - Computes the norm of the residual vector associated with
780 dsic.upv.es!jroman 935 
   the i-th computed left eigenvector (only available in two-sided eigensolvers).
936 
 
937 
   Collective on EPS
938 
 
939 
   Input Parameter:
940 
.  eps - the eigensolver context
941 
.  i   - the solution index
942 
 
943 
   Output Parameter:
944 
.  norm - the residual norm, computed as ||y'A-ky'B||_2 where k is the
945 
   eigenvalue and y is the left eigenvector.
946 
   If k=0 then the residual norm is computed as ||y'A||_2.
947 
 
948 
   Notes:
1267 slepc 949 
   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
780 dsic.upv.es!jroman 950 
   Eigenpairs are indexed according to the ordering criterion established
951 
   with EPSSetWhichEigenpairs().
952 
 
953 
   Level: beginner
954 
 
955 
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
956 
@*/
2331 jroman 957 
PetscErrorCode EPSComputeResidualNormLeft(EPS eps,PetscInt i,PetscReal *norm)
780 dsic.upv.es!jroman 958 
{
959 
  PetscErrorCode ierr;
2331 jroman 960 
  Vec            u,v,w,xr,xi;
961 
  Mat            A,B;
962 
  PetscScalar    kr,ki;
2320 jroman 963 
#if !defined(PETSC_USE_COMPLEX)
2331 jroman 964 
  PetscReal      ni,nr;
780 dsic.upv.es!jroman 965 
#endif
966 
 
967 
  PetscFunctionBegin;
2213 jroman 968 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2328 jroman 969 
  PetscValidLogicalCollectiveInt(eps,i,2);
970 
  PetscValidPointer(norm,3);
2762 jroman 971 
  if (!eps->leftvecs) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetLeftVectorsWanted");
780 dsic.upv.es!jroman 972 
  ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
2330 jroman 973 
  ierr = VecDuplicate(eps->W[0],&u);CHKERRQ(ierr);
974 
  ierr = VecDuplicate(eps->W[0],&v);CHKERRQ(ierr);
975 
  ierr = VecDuplicate(eps->W[0],&w);CHKERRQ(ierr);
976 
  ierr = VecDuplicate(eps->W[0],&xr);CHKERRQ(ierr);
977 
  ierr = VecDuplicate(eps->W[0],&xi);CHKERRQ(ierr);
978 
  ierr = EPSGetEigenvalue(eps,i,&kr,&ki);CHKERRQ(ierr);
979 
  ierr = EPSGetEigenvectorLeft(eps,i,xr,xi);CHKERRQ(ierr);
780 dsic.upv.es!jroman 980 
 
2320 jroman 981 
#if !defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 982 
  if (ki == 0 ||
983 
    PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
984 
#endif
2331 jroman 985 
    ierr = MatMultTranspose(A,xr,u);CHKERRQ(ierr); /* u=A'*x */
780 dsic.upv.es!jroman 986 
    if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
2331 jroman 987 
      if (eps->isgeneralized) { ierr = MatMultTranspose(B,xr,w);CHKERRQ(ierr); }
988 
      else { ierr = VecCopy(xr,w);CHKERRQ(ierr); } /* w=B'*x */
989 
      ierr = VecAXPY(u,-kr,w);CHKERRQ(ierr); /* u=A'*x-k*B'*x */
780 dsic.upv.es!jroman 990 
    }
2331 jroman 991 
    ierr = VecNorm(u,NORM_2,norm);CHKERRQ(ierr);  
2320 jroman 992 
#if !defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 993 
  } else {
2331 jroman 994 
    ierr = MatMultTranspose(A,xr,u);CHKERRQ(ierr); /* u=A'*xr */
995 
    if (eps->isgeneralized) { ierr = MatMultTranspose(B,xr,v);CHKERRQ(ierr); }
996 
    else { ierr = VecCopy(xr,v);CHKERRQ(ierr); } /* v=B'*xr */
997 
    ierr = VecAXPY(u,-kr,v);CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
998 
    if (eps->isgeneralized) { ierr = MatMultTranspose(B,xi,w);CHKERRQ(ierr); }
999 
    else { ierr = VecCopy(xi,w);CHKERRQ(ierr); } /* w=B'*xi */
1000 
    ierr = VecAXPY(u,ki,w);CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr+ki*B'*xi */
1001 
    ierr = VecNorm(u,NORM_2,&nr);CHKERRQ(ierr);
1002 
    ierr = MatMultTranspose(A,xi,u);CHKERRQ(ierr); /* u=A'*xi */
1003 
    ierr = VecAXPY(u,-kr,w);CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi */
1004 
    ierr = VecAXPY(u,-ki,v);CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi-ki*B'*xr */
1005 
    ierr = VecNorm(u,NORM_2,&ni);CHKERRQ(ierr);
1006 
    *norm = SlepcAbsEigenvalue(nr,ni);
780 dsic.upv.es!jroman 1007 
  }
1008 
#endif
1009 
 
2330 jroman 1010 
  ierr = VecDestroy(&w);CHKERRQ(ierr);
1011 
  ierr = VecDestroy(&v);CHKERRQ(ierr);
1012 
  ierr = VecDestroy(&u);CHKERRQ(ierr);
1013 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
1014 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1015 
  PetscFunctionReturn(0);
1016 
}
1017 
 
1018 
#undef __FUNCT__  
1812 antodo 1019 
#define __FUNCT__ "EPSComputeRelativeError_Private"
1020 
/*
1021 
   EPSComputeRelativeError_Private - Computes the relative error bound
1022 
   associated with an eigenpair.
1023 
*/
2331 jroman 1024 
PetscErrorCode EPSComputeRelativeError_Private(EPS eps,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,PetscReal *error)
528 dsic.upv.es!antodo 1025 
{
1026 
  PetscErrorCode ierr;
2331 jroman 1027 
  PetscReal      norm,er;
2320 jroman 1028 
#if !defined(PETSC_USE_COMPLEX)
528 dsic.upv.es!antodo 1029 
  PetscReal      ei;
1030 
#endif
1031 
 
1032 
  PetscFunctionBegin;
1893 jroman 1033 
  ierr = EPSComputeResidualNorm_Private(eps,kr,ki,xr,xi,&norm);CHKERRQ(ierr);
528 dsic.upv.es!antodo 1034 
 
2320 jroman 1035 
#if !defined(PETSC_USE_COMPLEX)
2115 eromero 1036 
  if (ki == 0) {
528 dsic.upv.es!antodo 1037 
#endif
1893 jroman 1038 
    ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
2320 jroman 1039 
#if !defined(PETSC_USE_COMPLEX)
528 dsic.upv.es!antodo 1040 
  } else {
1893 jroman 1041 
    ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);  
2115 eromero 1042 
    ierr = VecNorm(xi,NORM_2,&ei);CHKERRQ(ierr);
1043 
    er = SlepcAbsEigenvalue(er,ei);
528 dsic.upv.es!antodo 1044 
  }
1045 
#endif    
2219 jroman 1046 
  ierr = (*eps->conv_func)(eps,kr,ki,norm/er,error,eps->conv_ctx);CHKERRQ(ierr);
1812 antodo 1047 
  PetscFunctionReturn(0);
1048 
}
1049 
 
1050 
#undef __FUNCT__  
1051 
#define __FUNCT__ "EPSComputeRelativeError"
1052 
/*@
1053 
   EPSComputeRelativeError - Computes the relative error bound associated
1054 
   with the i-th computed eigenpair.
1055 
 
1056 
   Collective on EPS
1057 
 
1058 
   Input Parameter:
1059 
.  eps - the eigensolver context
1060 
.  i   - the solution index
1061 
 
1062 
   Output Parameter:
1063 
.  error - the relative error bound, computed as ||Ax-kBx||_2/||kx||_2 where
1064 
   k is the eigenvalue and x is the eigenvector.
1065 
   If k=0 the relative error is computed as ||Ax||_2/||x||_2.
1066 
 
1067 
   Level: beginner
1068 
 
1069 
.seealso: EPSSolve(), EPSComputeResidualNorm(), EPSGetErrorEstimate()
1070 
@*/
2331 jroman 1071 
PetscErrorCode EPSComputeRelativeError(EPS eps,PetscInt i,PetscReal *error)
1812 antodo 1072 
{
1073 
  PetscErrorCode ierr;
2331 jroman 1074 
  Vec            xr,xi;  
1075 
  PetscScalar    kr,ki;  
1812 antodo 1076 
 
1077 
  PetscFunctionBegin;
2213 jroman 1078 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);  
2328 jroman 1079 
  PetscValidLogicalCollectiveInt(eps,i,2);
1812 antodo 1080 
  PetscValidPointer(error,3);
2330 jroman 1081 
  ierr = VecDuplicate(eps->V[0],&xr);CHKERRQ(ierr);
1082 
  ierr = VecDuplicate(eps->V[0],&xi);CHKERRQ(ierr);
1083 
  ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
1084 
  ierr = EPSComputeRelativeError_Private(eps,kr,ki,xr,xi,error);CHKERRQ(ierr);  
1085 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
1086 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
528 dsic.upv.es!antodo 1087 
  PetscFunctionReturn(0);
1088 
}
1089 
 
1090 
#undef __FUNCT__  
780 dsic.upv.es!jroman 1091 
#define __FUNCT__ "EPSComputeRelativeErrorLeft"
1092 
/*@
1093 
   EPSComputeRelativeErrorLeft - Computes the relative error bound associated
1094 
   with the i-th computed eigenvalue and left eigenvector (only available in
1095 
   two-sided eigensolvers).
1096 
 
1097 
   Collective on EPS
1098 
 
1099 
   Input Parameter:
1100 
.  eps - the eigensolver context
1101 
.  i   - the solution index
1102 
 
1103 
   Output Parameter:
1104 
.  error - the relative error bound, computed as ||y'A-ky'B||_2/||ky||_2 where
1105 
   k is the eigenvalue and y is the left eigenvector.
1106 
   If k=0 the relative error is computed as ||y'A||_2/||y||_2.
1107 
 
1108 
   Level: beginner
1109 
 
1110 
.seealso: EPSSolve(), EPSComputeResidualNormLeft(), EPSGetErrorEstimateLeft()
1111 
@*/
2331 jroman 1112 
PetscErrorCode EPSComputeRelativeErrorLeft(EPS eps,PetscInt i,PetscReal *error)
780 dsic.upv.es!jroman 1113 
{
1114 
  PetscErrorCode ierr;
2331 jroman 1115 
  Vec            xr,xi;  
1116 
  PetscScalar    kr,ki;  
1117 
  PetscReal      norm,er;
2320 jroman 1118 
#if !defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 1119 
  Vec            u;
1120 
  PetscReal      ei;
1121 
#endif
1122 
 
1123 
  PetscFunctionBegin;
2213 jroman 1124 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);  
2330 jroman 1125 
  ierr = EPSComputeResidualNormLeft(eps,i,&norm);CHKERRQ(ierr);
2328 jroman 1126 
  PetscValidLogicalCollectiveInt(eps,i,2);
1127 
  PetscValidPointer(error,3);
2330 jroman 1128 
  ierr = VecDuplicate(eps->W[0],&xr);CHKERRQ(ierr);
1129 
  ierr = VecDuplicate(eps->W[0],&xi);CHKERRQ(ierr);
1130 
  ierr = EPSGetEigenvalue(eps,i,&kr,&ki);CHKERRQ(ierr);
1131 
  ierr = EPSGetEigenvectorLeft(eps,i,xr,xi);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1132 
 
2320 jroman 1133 
#if !defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 1134 
  if (ki == 0 ||
1135 
    PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
1136 
#endif
2331 jroman 1137 
    ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1138 
    if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
868 dsic.upv.es!antodo 1139 
      *error =  norm / (PetscAbsScalar(kr) * er);
1140 
    } else {
1141 
      *error = norm / er;
780 dsic.upv.es!jroman 1142 
    }
2320 jroman 1143 
#if !defined(PETSC_USE_COMPLEX)
780 dsic.upv.es!jroman 1144 
  } else {
2331 jroman 1145 
    ierr = VecDuplicate(xi,&u);CHKERRQ(ierr);  
1146 
    ierr = VecCopy(xi,u);CHKERRQ(ierr);  
1147 
    ierr = VecAXPBY(u,kr,-ki,xr);CHKERRQ(ierr);  
1148 
    ierr = VecNorm(u,NORM_2,&er);CHKERRQ(ierr);  
1149 
    ierr = VecAXPBY(xi,kr,ki,xr);CHKERRQ(ierr);      
1150 
    ierr = VecNorm(xi,NORM_2,&ei);CHKERRQ(ierr);  
2330 jroman 1151 
    ierr = VecDestroy(&u);CHKERRQ(ierr);  
2331 jroman 1152 
    *error = norm / SlepcAbsEigenvalue(er,ei);
780 dsic.upv.es!jroman 1153 
  }
1154 
#endif    
1155 
 
2330 jroman 1156 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
1157 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1158 
  PetscFunctionReturn(0);
1159 
}
1160 
 
677 dsic.upv.es!antodo 1161 
#define SWAP(a,b,t) {t=a;a=b;b=t;}
1162 
 
528 dsic.upv.es!antodo 1163 
#undef __FUNCT__  
677 dsic.upv.es!antodo 1164 
#define __FUNCT__ "EPSSortEigenvalues"
528 dsic.upv.es!antodo 1165 
/*@
1782 antodo 1166 
   EPSSortEigenvalues - Sorts a list of eigenvalues according to the criterion
1811 jroman 1167 
   specified via EPSSetWhichEigenpairs().
528 dsic.upv.es!antodo 1168 
 
677 dsic.upv.es!antodo 1169 
   Not Collective
528 dsic.upv.es!antodo 1170 
 
677 dsic.upv.es!antodo 1171 
   Input Parameters:
1783 antodo 1172 
+  eps   - the eigensolver context
1811 jroman 1173 
.  n     - number of eigenvalues in the list
1174 
.  eigr  - pointer to the array containing the eigenvalues
1782 antodo 1175 
-  eigi  - imaginary part of the eigenvalues (only when using real numbers)
528 dsic.upv.es!antodo 1176 
 
677 dsic.upv.es!antodo 1177 
   Output Parameter:
1811 jroman 1178 
.  perm  - resulting permutation
528 dsic.upv.es!antodo 1179 
 
1811 jroman 1180 
   Note:
677 dsic.upv.es!antodo 1181 
   The result is a list of indices in the original eigenvalue array
1182 
   corresponding to the first nev eigenvalues sorted in the specified
1811 jroman 1183 
   criterion.
528 dsic.upv.es!antodo 1184 
 
677 dsic.upv.es!antodo 1185 
   Level: developer
528 dsic.upv.es!antodo 1186 
 
1628 slepc 1187 
.seealso: EPSSortEigenvaluesReal(), EPSSetWhichEigenpairs()
528 dsic.upv.es!antodo 1188 
@*/
1782 antodo 1189 
PetscErrorCode EPSSortEigenvalues(EPS eps,PetscInt n,PetscScalar *eigr,PetscScalar *eigi,PetscInt *perm)
528 dsic.upv.es!antodo 1190 
{
1191 
  PetscErrorCode ierr;
1782 antodo 1192 
  PetscScalar    re,im;
1193 
  PetscInt       i,j,result,tmp;
528 dsic.upv.es!antodo 1194 
 
1195 
  PetscFunctionBegin;
2349 jroman 1196 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);  
1197 
  PetscValidScalarPointer(eigr,3);
1198 
  PetscValidScalarPointer(eigi,4);
1199 
  PetscValidIntPointer(perm,5);
1782 antodo 1200 
  for (i=0; i<n; i++) { perm[i] = i; }
1201 
  /* insertion sort */
1783 antodo 1202 
  for (i=n-1; i>=0; i--) {
1782 antodo 1203 
    re = eigr[perm[i]];
1834 antodo 1204 
    im = eigi[perm[i]];
1783 antodo 1205 
    j = i + 1;
2320 jroman 1206 
#if !defined(PETSC_USE_COMPLEX)
1783 antodo 1207 
    if (im != 0) {
1208 
      /* complex eigenvalue */
1209 
      i--;
1210 
      im = eigi[perm[i]];
1211 
    }
1212 
#endif
1213 
    while (j<n) {
1214 
      ierr = EPSCompareEigenvalues(eps,re,im,eigr[perm[j]],eigi[perm[j]],&result);CHKERRQ(ierr);
2097 eromero 1215 
      if (result < 0) break;
2320 jroman 1216 
#if !defined(PETSC_USE_COMPLEX)
1782 antodo 1217 
      /* keep together every complex conjugated eigenpair */
1218 
      if (im == 0) {
1219 
        if (eigi[perm[j]] == 0) {
1783 antodo 1220 
#endif
1221 
          tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
1222 
          j++;
2320 jroman 1223 
#if !defined(PETSC_USE_COMPLEX)
1782 antodo 1224 
        } else {
1783 antodo 1225 
          tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
1226 
          j+=2;
1782 antodo 1227 
        }
1783 antodo 1228 
      } else {
1782 antodo 1229 
        if (eigi[perm[j]] == 0) {
1783 antodo 1230 
          tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = perm[j-1]; perm[j-1] = tmp;
1231 
          j++;
1782 antodo 1232 
        } else {
1783 antodo 1233 
          tmp = perm[j-2]; perm[j-2] = perm[j]; perm[j] = tmp;
1782 antodo 1234 
          tmp = perm[j-1]; perm[j-1] = perm[j+1]; perm[j+1] = tmp;
1783 antodo 1235 
          j+=2;
1782 antodo 1236 
        }
1237 
      }
677 dsic.upv.es!antodo 1238 
#endif
1239 
    }
1240 
  }
528 dsic.upv.es!antodo 1241 
  PetscFunctionReturn(0);
1242 
}
689 dsic.upv.es!jroman 1243 
 
1244 
#undef __FUNCT__  
1477 slepc 1245 
#define __FUNCT__ "EPSSortEigenvaluesReal"
1246 
/*@
1247 
   EPSSortEigenvaluesReal - Sorts a list of eigenvalues according to a certain
1248 
   criterion (version for real eigenvalues only).
1249 
 
1250 
   Not Collective
1251 
 
1252 
   Input Parameters:
1811 jroman 1253 
+  eps   - the eigensolver context
1254 
.  n     - number of eigenvalue in the list
1255 
-  eig   - pointer to the array containing the eigenvalues (real)
1477 slepc 1256 
 
1257 
   Output Parameter:
1811 jroman 1258 
.  perm  - resulting permutation
1477 slepc 1259 
 
1811 jroman 1260 
   Note:
1477 slepc 1261 
   The result is a list of indices in the original eigenvalue array
1262 
   corresponding to the first nev eigenvalues sorted in the specified
1811 jroman 1263 
   criterion.
1477 slepc 1264 
 
1265 
   Level: developer
1266 
 
1811 jroman 1267 
.seealso: EPSSortEigenvalues(), EPSSetWhichEigenpairs(), EPSCompareEigenvalues()
1477 slepc 1268 
@*/
1782 antodo 1269 
PetscErrorCode EPSSortEigenvaluesReal(EPS eps,PetscInt n,PetscReal *eig,PetscInt *perm)
1477 slepc 1270 
{
1271 
  PetscErrorCode ierr;
1782 antodo 1272 
  PetscScalar    re;
1273 
  PetscInt       i,j,result,tmp;
1477 slepc 1274 
 
1275 
  PetscFunctionBegin;
2349 jroman 1276 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);  
1277 
  PetscValidPointer(eig,3);
1278 
  PetscValidIntPointer(perm,4);
1782 antodo 1279 
  for (i=0; i<n; i++) { perm[i] = i; }
1280 
  /* insertion sort */
1281 
  for (i=1; i<n; i++) {
1282 
    re = eig[perm[i]];
1283 
    j = i-1;
1284 
    ierr = EPSCompareEigenvalues(eps,re,0.0,eig[perm[j]],0.0,&result);CHKERRQ(ierr);
2097 eromero 1285 
    while (result<=0 && j>=0) {
1782 antodo 1286 
      tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
1287 
      if (j>=0) {
1288 
        ierr = EPSCompareEigenvalues(eps,re,0.0,eig[perm[j]],0.0,&result);CHKERRQ(ierr);
1289 
      }
1290 
    }
1291 
  }
1292 
  PetscFunctionReturn(0);
1293 
}
1477 slepc 1294 
 
1782 antodo 1295 
#undef __FUNCT__  
1296 
#define __FUNCT__ "EPSCompareEigenvalues"
1811 jroman 1297 
/*@
1298 
   EPSCompareEigenvalues - Compares two (possibly complex) eigenvalues according
1299 
   to a certain criterion.
1300 
 
1301 
   Not Collective
1302 
 
1303 
   Input Parameters:
1304 
+  eps   - the eigensolver context
1305 
.  ar     - real part of the 1st eigenvalue
1306 
.  ai     - imaginary part of the 1st eigenvalue
1307 
.  br     - real part of the 2nd eigenvalue
1308 
-  bi     - imaginary part of the 2nd eigenvalue
1309 
 
1310 
   Output Parameter:
1311 
.  res    - result of comparison
1312 
 
1313 
   Notes:
2097 eromero 1314 
   The returning parameter 'res' can be:
1315 
+  negative - if the 1st eigenvalue is preferred to the 2st one
1316 
.  zero     - if both eigenvalues are equally preferred
1317 
-  positive - if the 2st eigenvalue is preferred to the 1st one
1811 jroman 1318 
 
1319 
   The criterion of comparison is related to the 'which' parameter set with
1320 
   EPSSetWhichEigenpairs().
1321 
 
1322 
   Level: developer
1323 
 
1324 
.seealso: EPSSortEigenvalues(), EPSSetWhichEigenpairs()
1325 
@*/
1782 antodo 1326 
PetscErrorCode EPSCompareEigenvalues(EPS eps,PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *result)
1327 
{
1328 
  PetscErrorCode ierr;
1329 
 
1330 
  PetscFunctionBegin;
2349 jroman 1331 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);  
1332 
  PetscValidIntPointer(result,6);
2760 jroman 1333 
  if (!eps->which_func) SETERRQ(((PetscObject)eps)->comm,1,"Undefined eigenvalue comparison function");
1334 
  ierr = (*eps->which_func)(ar,ai,br,bi,result,eps->which_ctx);CHKERRQ(ierr);
1477 slepc 1335 
  PetscFunctionReturn(0);
1336 
}
1337 
 
1338 
#undef __FUNCT__  
689 dsic.upv.es!jroman 1339 
#define __FUNCT__ "EPSGetStartVector"
1340 
/*@
1937 jroman 1341 
   EPSGetStartVector - Gets a suitable vector to be used as the starting vector
1342 
   for the recurrence that builds the right subspace.
689 dsic.upv.es!jroman 1343 
 
1344 
   Collective on EPS and Vec
1345 
 
1346 
   Input Parameters:
1347 
+  eps - the eigensolver context
1937 jroman 1348 
-  i   - iteration number
689 dsic.upv.es!jroman 1349 
 
1059 slepc 1350 
   Output Parameters:
1351 
+  vec - the start vector
1352 
-  breakdown - flag indicating that a breakdown has occurred
689 dsic.upv.es!jroman 1353 
 
1354 
   Notes:
1355 
   The start vector is computed from another vector: for the first step (i=0),
1933 jroman 1356 
   the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
1229 slepc 1357 
   vector is created. Then this vector is forced to be in the range of OP (only
1358 
   for generalized definite problems) and orthonormalized with respect to all
1359 
   V-vectors up to i-1.
689 dsic.upv.es!jroman 1360 
 
1059 slepc 1361 
   The flag breakdown is set to true if either i=0 and the vector belongs to the
1362 
   deflation space, or i>0 and the vector is linearly dependent with respect
1363 
   to the V-vectors.
1364 
 
689 dsic.upv.es!jroman 1365 
   The caller must pass a vector already allocated with dimensions conforming
1366 
   to the initial vector. This vector is overwritten.
1367 
 
1368 
   Level: developer
1369 
 
1933 jroman 1370 
.seealso: EPSSetInitialSpace()
689 dsic.upv.es!jroman 1371 
@*/
2216 jroman 1372 
PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,Vec vec,PetscBool *breakdown)
689 dsic.upv.es!jroman 1373 
{
1374 
  PetscErrorCode ierr;
1375 
  PetscReal      norm;
2216 jroman 1376 
  PetscBool      lindep;
689 dsic.upv.es!jroman 1377 
  Vec            w;
1378 
 
1379 
  PetscFunctionBegin;
2213 jroman 1380 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 1381 
  PetscValidLogicalCollectiveInt(eps,i,2);
2213 jroman 1382 
  PetscValidHeaderSpecific(vec,VEC_CLASSID,3);
2326 jroman 1383 
  PetscCheckSameComm(eps,1,vec,3);
689 dsic.upv.es!jroman 1384 
 
1933 jroman 1385 
  ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
1386 
 
1387 
  /* For the first step, use the first initial vector, otherwise a random one */
1388 
  if (i==0 && eps->nini>0) {
1389 
    ierr = VecCopy(eps->V[0],w);CHKERRQ(ierr);
1057 slepc 1390 
  } else {
2027 jroman 1391 
    ierr = SlepcVecSetRandom(w,eps->rand);CHKERRQ(ierr);
689 dsic.upv.es!jroman 1392 
  }
1393 
 
1229 slepc 1394 
  /* Force the vector to be in the range of OP for definite generalized problems */
2743 carcamgo 1395 
  if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
1229 slepc 1396 
    ierr = STApply(eps->OP,w,vec);CHKERRQ(ierr);
1397 
  } else {
1398 
    ierr = VecCopy(w,vec);CHKERRQ(ierr);
1399 
  }
689 dsic.upv.es!jroman 1400 
 
1401 
  /* Orthonormalize the vector with respect to previous vectors */
1755 antodo 1402 
  ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,i,PETSC_NULL,eps->V,vec,PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
1057 slepc 1403 
  if (breakdown) *breakdown = lindep;
1169 slepc 1404 
  else if (lindep || norm == 0.0) {
2762 jroman 1405 
    if (i==0) SETERRQ(((PetscObject)eps)->comm,1,"Initial vector is zero or belongs to the deflation space");
1406 
    else SETERRQ(((PetscObject)eps)->comm,1,"Unable to generate more start vectors");
750 dsic.upv.es!antodo 1407 
  }
1509 slepc 1408 
  ierr = VecScale(vec,1.0/norm);CHKERRQ(ierr);
689 dsic.upv.es!jroman 1409 
 
2305 jroman 1410 
  ierr = VecDestroy(&w);CHKERRQ(ierr);
689 dsic.upv.es!jroman 1411 
  PetscFunctionReturn(0);
1412 
}
1936 jroman 1413 
 
780 dsic.upv.es!jroman 1414 
#undef __FUNCT__  
1936 jroman 1415 
#define __FUNCT__ "EPSGetStartVectorLeft"
780 dsic.upv.es!jroman 1416 
/*@
1937 jroman 1417 
   EPSGetStartVectorLeft - Gets a suitable vector to be used as the starting vector
1418 
   in the recurrence that builds the left subspace (in methods that work with two
1419 
   subspaces).
689 dsic.upv.es!jroman 1420 
 
780 dsic.upv.es!jroman 1421 
   Collective on EPS and Vec
1422 
 
1423 
   Input Parameters:
1424 
+  eps - the eigensolver context
1937 jroman 1425 
-  i   - iteration number
780 dsic.upv.es!jroman 1426 
 
1427 
   Output Parameter:
1937 jroman 1428 
+  vec - the start vector
1429 
-  breakdown - flag indicating that a breakdown has occurred
780 dsic.upv.es!jroman 1430 
 
1431 
   Notes:
1432 
   The start vector is computed from another vector: for the first step (i=0),
1937 jroman 1433 
   the first left initial vector is used (see EPSSetInitialSpaceLeft()); otherwise
780 dsic.upv.es!jroman 1434 
   a random vector is created. Then this vector is forced to be in the range
1435 
   of OP' and orthonormalized with respect to all W-vectors up to i-1.
1436 
 
1937 jroman 1437 
   The flag breakdown is set to true if i>0 and the vector is linearly dependent
1438 
   with respect to the W-vectors.
1439 
 
780 dsic.upv.es!jroman 1440 
   The caller must pass a vector already allocated with dimensions conforming
1441 
   to the left initial vector. This vector is overwritten.
1442 
 
1443 
   Level: developer
1444 
 
1937 jroman 1445 
.seealso: EPSSetInitialSpaceLeft()
780 dsic.upv.es!jroman 1446 
 
1447 
@*/
2216 jroman 1448 
PetscErrorCode EPSGetStartVectorLeft(EPS eps,PetscInt i,Vec vec,PetscBool *breakdown)
780 dsic.upv.es!jroman 1449 
{
1450 
  PetscErrorCode ierr;
1451 
  PetscReal      norm;
2216 jroman 1452 
  PetscBool      lindep;
780 dsic.upv.es!jroman 1453 
  Vec            w;
1454 
 
1455 
  PetscFunctionBegin;
2213 jroman 1456 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 1457 
  PetscValidLogicalCollectiveInt(eps,i,2);
2213 jroman 1458 
  PetscValidHeaderSpecific(vec,VEC_CLASSID,3);
2326 jroman 1459 
  PetscCheckSameComm(eps,1,vec,3);
780 dsic.upv.es!jroman 1460 
 
1937 jroman 1461 
  ierr = VecDuplicate(eps->W[0],&w);CHKERRQ(ierr);
1462 
 
1463 
  /* For the first step, use the first initial left vector, otherwise a random one */
1464 
  if (i==0 && eps->ninil>0) {
1465 
    ierr = VecCopy(eps->W[0],w);CHKERRQ(ierr);
1466 
  } else {
2027 jroman 1467 
    ierr = SlepcVecSetRandom(w,eps->rand);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1468 
  }
1469 
 
1937 jroman 1470 
  /* Force the vector to be in the range of OP' */
780 dsic.upv.es!jroman 1471 
  ierr = STApplyTranspose(eps->OP,w,vec);CHKERRQ(ierr);
1472 
 
1473 
  /* Orthonormalize the vector with respect to previous vectors */
1937 jroman 1474 
  ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,i,PETSC_NULL,eps->W,vec,PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
1475 
  if (breakdown) *breakdown = lindep;
1476 
  else if (lindep || norm == 0.0) {
2762 jroman 1477 
    if (i==0) SETERRQ(((PetscObject)eps)->comm,1,"Left initial vector is zero");
1478 
    else SETERRQ(((PetscObject)eps)->comm,1,"Unable to generate more left start vectors");
780 dsic.upv.es!jroman 1479 
  }
828 dsic.upv.es!antodo 1480 
  ierr = VecScale(vec,1/norm);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1481 
 
2305 jroman 1482 
  ierr = VecDestroy(&w);CHKERRQ(ierr);
780 dsic.upv.es!jroman 1483 
  PetscFunctionReturn(0);
1484 
}