Subversion Repositories slepc-dev

Rev

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