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