Subversion Repositories slepc-dev

Rev

Details | Compare with Previous | Last modification | View Log | RSS feed

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