Subversion Repositories slepc-dev

Rev

Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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