Subversion Repositories slepc-dev

/*

      EPS routines related to the solution process.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

#include "private/epsimpl.h"   /*I "slepceps.h" I*/

#undef __FUNCT__  

#define __FUNCT__ "EPSSolve"

/*@

   EPSSolve - Solves the eigensystem.

   Collective on EPS

   Input Parameter:

.  eps - eigensolver context obtained from EPSCreate()

   Options Database:

+   -eps_view - print information about the solver used

.   -eps_view_binary - save the matrices to the default binary file

-   -eps_plot_eigs - plot computed eigenvalues

   Level: beginner

.seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()

@*/

PetscErrorCode EPSSolve(EPS eps)

{

  PetscErrorCode ierr;

  PetscInt       i;

  PetscReal      re,im;

  PetscTruth     flg;

  PetscViewer    viewer;

  PetscDraw      draw;

  PetscDrawSP    drawsp;

  STMatMode      matmode;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_view_binary",&flg);CHKERRQ(ierr);

  if (flg) {

    Mat A,B;

    ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);

    ierr = MatView(A,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));CHKERRQ(ierr);

    if (B) ierr = MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));CHKERRQ(ierr);

  }

  /* reset the convergence flag from the previous solves */

  eps->reason = EPS_CONVERGED_ITERATING;

  if (!eps->setupcalled){ ierr = EPSSetUp(eps);CHKERRQ(ierr); }

  ierr = STResetOperationCounters(eps->OP);CHKERRQ(ierr);

  ierr = IPResetOperationCounters(eps->ip);CHKERRQ(ierr);

  eps->evecsavailable = PETSC_FALSE;

  eps->nconv = 0;

  eps->its = 0;

  for (i=0;i<eps->ncv;i++) eps->eigr[i]=eps->eigi[i]=eps->errest[i]=0.0;

  EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->ncv);

  ierr = PetscLogEventBegin(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);

  switch (eps->solverclass) {

    case EPS_ONE_SIDE:

      ierr = (*eps->ops->solve)(eps);CHKERRQ(ierr); break;

    case EPS_TWO_SIDE:

      if (eps->ops->solvets) {

        ierr = (*eps->ops->solvets)(eps);CHKERRQ(ierr); break;

      } else {

        SETERRQ(1,"Two-sided version unavailable for this solver");

      }

    default:

      SETERRQ(1,"Wrong value of eps->solverclass");

  }

  ierr = STGetMatMode(eps->OP,&matmode);CHKERRQ(ierr);

  if (matmode == STMATMODE_INPLACE && eps->ispositive) {

    /* Purge eigenvectors before reverting operator */

    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);    

  }

  ierr = STPostSolve(eps->OP);CHKERRQ(ierr);

  ierr = PetscLogEventEnd(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);

  if (!eps->reason) {

    SETERRQ(1,"Internal error, solver returned without setting converged reason");

  }

  /* Map eigenvalues back to the original problem, necessary in some

  * spectral transformations */

  ierr = (*eps->ops->backtransform)(eps);CHKERRQ(ierr);

  /* Adjust left eigenvectors in generalized problems: y = B^T y */

  if (eps->isgeneralized && eps->solverclass == EPS_TWO_SIDE) {

    Mat B;

    KSP ksp;

    Vec w;

    ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);

    ierr = KSPCreate(((PetscObject)eps)->comm,&ksp);CHKERRQ(ierr);

    ierr = KSPSetOperators(ksp,B,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

    ierr = MatGetVecs(B,PETSC_NULL,&w);CHKERRQ(ierr);

    for (i=0;i<eps->nconv;i++) {

      ierr = VecCopy(eps->W[i],w);CHKERRQ(ierr);

      ierr = KSPSolveTranspose(ksp,w,eps->W[i]);CHKERRQ(ierr);

    }

    ierr = KSPDestroy(ksp);CHKERRQ(ierr);

    ierr = VecDestroy(w);CHKERRQ(ierr);

  }

#ifndef PETSC_USE_COMPLEX

  /* reorder conjugate eigenvalues (positive imaginary first) */

  for (i=0; i<eps->nconv-1; i++) {

    if (eps->eigi[i] != 0) {

      if (eps->eigi[i] < 0) {

        eps->eigi[i] = -eps->eigi[i];

        eps->eigi[i+1] = -eps->eigi[i+1];

        ierr = VecScale(eps->V[i+1],-1.0); CHKERRQ(ierr);

      }

      i++;

    }

  }

#endif

  /* sort eigenvalues according to eps->which parameter */

  ierr = PetscFree(eps->perm);CHKERRQ(ierr);

  if (eps->nconv > 0) {

    ierr = PetscMalloc(sizeof(PetscInt)*eps->nconv, &eps->perm); CHKERRQ(ierr);

    ierr = EPSSortEigenvalues(eps->nconv, eps->eigr, eps->eigi, eps->which, eps->nconv, eps->perm); CHKERRQ(ierr);

  }

  ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_view",&flg);CHKERRQ(ierr);

  if (flg && !PetscPreLoadingOn) { ierr = EPSView(eps,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); }

  ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_plot_eigs",&flg);CHKERRQ(ierr);

  if (flg) {

    ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Computed Eigenvalues",

                             PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);CHKERRQ(ierr);

    ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);

    ierr = PetscDrawSPCreate(draw,1,&drawsp);CHKERRQ(ierr);

    for( i=0; i<eps->nconv; i++ ) {

#if defined(PETSC_USE_COMPLEX)

      re = PetscRealPart(eps->eigr[i]);

      im = PetscImaginaryPart(eps->eigi[i]);

#else

      re = eps->eigr[i];

      im = eps->eigi[i];

#endif

      ierr = PetscDrawSPAddPoint(drawsp,&re,&im);CHKERRQ(ierr);

    }

    ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr);

    ierr = PetscDrawSPDestroy(drawsp);CHKERRQ(ierr);

    ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetIterationNumber"

/*@

   EPSGetIterationNumber - Gets the current iteration number. If the

   call to EPSSolve() is complete, then it returns the number of iterations

   carried out by the solution method.

   Not Collective

   Input Parameter:

.  eps - the eigensolver context

   Output Parameter:

.  its - number of iterations

   Level: intermediate

   Note:

   During the i-th iteration this call returns i-1. If EPSSolve() is

   complete, then parameter "its" contains either the iteration number at

   which convergence was successfully reached, or failure was detected.  

   Call EPSGetConvergedReason() to determine if the solver converged or

   failed and why.

.seealso: EPSGetConvergedReason(), EPSSetTolerances()

@*/

PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  PetscValidIntPointer(its,2);

  *its = eps->its;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetOperationCounters"

/*@

   EPSGetOperationCounters - Gets the total number of operator applications,

   inner product operations and linear iterations used by the ST object

   during the last EPSSolve() call.

   Not Collective

   Input Parameter:

.  eps - EPS context

   Output Parameter:

+  ops  - number of operator applications

.  dots - number of inner product operations

-  lits - number of linear iterations

   Notes:

   When the eigensolver algorithm invokes STApply() then a linear system

   must be solved (except in the case of standard eigenproblems and shift

   transformation). The number of iterations required in this solve is

   accumulated into a counter whose value is returned by this function.

   These counters are reset to zero at each successive call to EPSSolve().

   Level: intermediate

@*/

PetscErrorCode EPSGetOperationCounters(EPS eps,PetscInt* ops,PetscInt* dots,PetscInt* lits)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  ierr = STGetOperationCounters(eps->OP,ops,lits);CHKERRQ(ierr);

  if (dots) {

    ierr = IPGetOperationCounters(eps->ip,dots);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetConverged"

/*@

   EPSGetConverged - Gets the number of converged eigenpairs.

   Not Collective

   Input Parameter:

.  eps - the eigensolver context

   Output Parameter:

.  nconv - number of converged eigenpairs

   Note:

   This function should be called after EPSSolve() has finished.

   Level: beginner

.seealso: EPSSetDimensions()

@*/

PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  PetscValidIntPointer(nconv,2);

  *nconv = eps->nconv;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetConvergedReason"

/*@C

   EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was

   stopped.

   Not Collective

   Input Parameter:

.  eps - the eigensolver context

   Output Parameter:

.  reason - negative value indicates diverged, positive value converged

   (see EPSConvergedReason)

   Possible values for reason:

+  EPS_CONVERGED_TOL - converged up to tolerance

.  EPS_DIVERGED_ITS - required more than its to reach convergence

.  EPS_DIVERGED_BREAKDOWN - generic breakdown in method

-  EPS_DIVERGED_NONSYMMETRIC - The operator is nonsymmetric

   Level: intermediate

   Notes: Can only be called after the call to EPSSolve() is complete.

.seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason

@*/

PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  PetscValidIntPointer(reason,2);

  *reason = eps->reason;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetInvariantSubspace" 

/*@

   EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant

   subspace.

   Not Collective

   Input Parameter:

.  eps - the eigensolver context

   Output Parameter:

.  v - an array of vectors

   Notes:

   This function should be called after EPSSolve() has finished.

   The user should provide in v an array of nconv vectors, where nconv is

   the value returned by EPSGetConverged().

   The first k vectors returned in v span an invariant subspace associated

   with the first k computed eigenvalues (note that this is not true if the

   k-th eigenvalue is complex and matrix A is real; in this case the first

   k+1 vectors should be used). An invariant subspace X of A satisfies Ax

   in X for all x in X (a similar definition applies for generalized

   eigenproblems).

   Level: intermediate

.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetLeftInvariantSubspace()

@*/

PetscErrorCode EPSGetInvariantSubspace(EPS eps, Vec *v)

{

  PetscErrorCode ierr;

  PetscInt       i;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  PetscValidPointer(v,2);

  PetscValidHeaderSpecific(*v,VEC_COOKIE,2);

  if (!eps->V) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

  }

  if (!eps->ishermitian && eps->evecsavailable) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetRightVector,EPSComputeRelativeError or EPSComputeResidualNorm");

  }

  for (i=0;i<eps->nconv;i++) {

    ierr = VecCopy(eps->V[i],v[i]);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetLeftInvariantSubspace" 

/*@

   EPSGetLeftInvariantSubspace - Gets an orthonormal basis of the computed left

   invariant subspace (only available in two-sided eigensolvers).

   Not Collective

   Input Parameter:

.  eps - the eigensolver context

   Output Parameter:

.  v - an array of vectors

   Notes:

   This function should be called after EPSSolve() has finished.

   The user should provide in v an array of nconv vectors, where nconv is

   the value returned by EPSGetConverged().

   The first k vectors returned in v span a left invariant subspace associated

   with the first k computed eigenvalues (note that this is not true if the

   k-th eigenvalue is complex and matrix A is real; in this case the first

   k+1 vectors should be used). A left invariant subspace Y of A satisfies y'A

   in Y for all y in Y (a similar definition applies for generalized

   eigenproblems).

   Level: intermediate

.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetInvariantSubspace

@*/

PetscErrorCode EPSGetLeftInvariantSubspace(EPS eps, Vec *v)

{

  PetscErrorCode ierr;

  PetscInt       i;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  PetscValidPointer(v,2);

  PetscValidHeaderSpecific(*v,VEC_COOKIE,2);

  if (!eps->W) {

    if (eps->solverclass!=EPS_TWO_SIDE) {

      SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");

    } else {

      SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

    }

  }

  for (i=0;i<eps->nconv;i++) {

    ierr = VecCopy(eps->W[i],v[i]);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetEigenpair" 

/*@

   EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by

   EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.

   Not Collective

   Input Parameters:

+  eps - eigensolver context

-  i   - index of the solution

   Output Parameters:

+  eigr - real part of eigenvalue

.  eigi - imaginary part of eigenvalue

.  Vr   - real part of eigenvector

-  Vi   - imaginary part of eigenvector

   Notes:

   If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is

   configured with complex scalars the eigenvalue is stored

   directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is

   set to zero).

   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).

   Eigenpairs are indexed according to the ordering criterion established

   with EPSSetWhichEigenpairs().

   Level: beginner

.seealso: EPSGetValue(), EPSGetRightVector(), EPSGetLeftVector(), EPSSolve(),

          EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()

@*/

PetscErrorCode EPSGetEigenpair(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi, Vec Vr, Vec Vi)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  if (!eps->eigr || !eps->eigi || !eps->V) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

  }

  if (i<0 || i>=eps->nconv) {

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");

  }

  ierr = EPSGetValue(eps,i,eigr,eigi);CHKERRQ(ierr);

  ierr = EPSGetRightVector(eps,i,Vr,Vi);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetValue" 

/*@

   EPSGetValue - Gets the i-th eigenvalue as computed by EPSSolve().

   Not Collective

   Input Parameters:

+  eps - eigensolver context

-  i   - index of the solution

   Output Parameters:

+  eigr - real part of eigenvalue

-  eigi - imaginary part of eigenvalue

   Notes:

   If the eigenvalue is real, then eigi is set to zero. If PETSc is

   configured with complex scalars the eigenvalue is stored

   directly in eigr (eigi is set to zero).

   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).

   Eigenpairs are indexed according to the ordering criterion established

   with EPSSetWhichEigenpairs().

   Level: beginner

.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),

          EPSGetEigenpair()

@*/

PetscErrorCode EPSGetValue(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi)

{

  PetscInt       k;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  if (!eps->eigr || !eps->eigi) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

  }

  if (i<0 || i>=eps->nconv) {

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");

  }

  if (!eps->perm) k = i;

  else k = eps->perm[i];

#ifdef PETSC_USE_COMPLEX

  if (eigr) *eigr = eps->eigr[k];

  if (eigi) *eigi = 0;

#else

  if (eigr) *eigr = eps->eigr[k];

  if (eigi) *eigi = eps->eigi[k];

#endif

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetRightVector" 

/*@

   EPSGetRightVector - Gets the i-th right eigenvector as computed by EPSSolve().

   Not Collective

   Input Parameters:

+  eps - eigensolver context

-  i   - index of the solution

   Output Parameters:

+  Vr   - real part of eigenvector

-  Vi   - imaginary part of eigenvector

   Notes:

   If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is

   configured with complex scalars the eigenvector is stored

   directly in Vr (Vi is set to zero).

   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).

   Eigenpairs are indexed according to the ordering criterion established

   with EPSSetWhichEigenpairs().

   Level: beginner

.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),

          EPSGetEigenpair(), EPSGetLeftVector()

@*/

PetscErrorCode EPSGetRightVector(EPS eps, PetscInt i, Vec Vr, Vec Vi)

{

  PetscErrorCode ierr;

  PetscInt       k;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  if (!eps->V) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

  }

  if (i<0 || i>=eps->nconv) {

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");

  }

  if (!eps->evecsavailable && (Vr || Vi) ) {

    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);

  }  

  if (!eps->perm) k = i;

  else k = eps->perm[i];

#ifdef PETSC_USE_COMPLEX

  if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }

  if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }

#else

  if (eps->eigi[k] > 0) { /* first value of conjugate pair */

    if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }

    if (Vi) { ierr = VecCopy(eps->V[k+1], Vi); CHKERRQ(ierr); }

  } else if (eps->eigi[k] < 0) { /* second value of conjugate pair */

    if (Vr) { ierr = VecCopy(eps->V[k-1], Vr); CHKERRQ(ierr); }

    if (Vi) {

      ierr = VecCopy(eps->V[k], Vi); CHKERRQ(ierr);

      ierr = VecScale(Vi,-1.0); CHKERRQ(ierr);

    }

  } else { /* real eigenvalue */

    if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }

    if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }

  }

#endif

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetLeftVector" 

/*@

   EPSGetLeftVector - Gets the i-th left eigenvector as computed by EPSSolve()

   (only available in two-sided eigensolvers).

   Not Collective

   Input Parameters:

+  eps - eigensolver context

-  i   - index of the solution

   Output Parameters:

+  Wr   - real part of eigenvector

-  Wi   - imaginary part of eigenvector

   Notes:

   If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is

   configured with complex scalars the eigenvector is stored

   directly in Wr (Wi is set to zero).

   The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).

   Eigenpairs are indexed according to the ordering criterion established

   with EPSSetWhichEigenpairs().

   Level: beginner

.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),

          EPSGetEigenpair(), EPSGetLeftVector()

@*/

PetscErrorCode EPSGetLeftVector(EPS eps, PetscInt i, Vec Wr, Vec Wi)

{

  PetscErrorCode ierr;

  PetscInt       k;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  if (!eps->W) {

    if (eps->solverclass!=EPS_TWO_SIDE) {

      SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");

    } else {

      SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

    }

  }

  if (i<0 || i>=eps->nconv) {

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");

  }

  if (!eps->evecsavailable && (Wr || Wi) ) {

    ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);

  }  

  if (!eps->perm) k = i;

  else k = eps->perm[i];

#ifdef PETSC_USE_COMPLEX

  if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }

  if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }

#else

  if (eps->eigi[k] > 0) { /* first value of conjugate pair */

    if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }

    if (Wi) { ierr = VecCopy(eps->W[k+1], Wi); CHKERRQ(ierr); }

  } else if (eps->eigi[k] < 0) { /* second value of conjugate pair */

    if (Wr) { ierr = VecCopy(eps->W[k-1], Wr); CHKERRQ(ierr); }

    if (Wi) {

      ierr = VecCopy(eps->W[k], Wi); CHKERRQ(ierr);

      ierr = VecScale(Wi,-1.0); CHKERRQ(ierr);

    }

  } else { /* real eigenvalue */

    if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }

    if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }

  }

#endif

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetErrorEstimate" 

/*@

   EPSGetErrorEstimate - Returns the error estimate associated to the i-th

   computed eigenpair.

   Not Collective

   Input Parameter:

+  eps - eigensolver context

-  i   - index of eigenpair

   Output Parameter:

.  errest - the error estimate

   Notes:

   This is the error estimate used internally by the eigensolver. The actual

   error bound can be computed with EPSComputeRelativeError(). See also the user's

   manual for details.

   Level: advanced

.seealso: EPSComputeRelativeError()

@*/

PetscErrorCode EPSGetErrorEstimate(EPS eps, PetscInt i, PetscReal *errest)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);

  if (!eps->eigr || !eps->eigi) {

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");

  }

  if (i<0 || i>=eps->nconv) {

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");

  }

  if (eps->perm) i = eps->perm[i];  

  if (errest) *errest = eps->errest[i];

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSGetErrorEstimateLeft" 

/*@

   EPSGetErrorEstimateLeft - Returns the left error estimate associated to the i-th

   computed eigenpair (only available in two-sided eigensolvers).

   Not Collective

   Input Parameter:

+  eps - eigensolver context

-  i   - index of eigenpair

   Output Parameter:

.  errest - the left error estimate

   Notes: