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