/*
|
/*
|
EPS routines related to the solution process.
|
EPS routines related to the solution process.
|
|
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
|
Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain
|
|
|
This file is part of SLEPc. See the README file for conditions of use
|
This file is part of SLEPc.
|
and additional information.
|
|
|
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*/
|
#include "private/epsimpl.h" /*I "slepceps.h" I*/
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSSolve"
|
#define __FUNCT__ "EPSSolve"
|
/*@
|
/*@
|
EPSSolve - Solves the eigensystem.
|
EPSSolve - Solves the eigensystem.
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - eigensolver context obtained from EPSCreate()
|
. eps - eigensolver context obtained from EPSCreate()
|
|
|
Options Database:
|
Options Database:
|
+ -eps_view - print information about the solver used
|
+ -eps_view - print information about the solver used
|
. -eps_view_binary - save the matrices to the default binary file
|
. -eps_view_binary - save the matrices to the default binary file
|
- -eps_plot_eigs - plot computed eigenvalues
|
- -eps_plot_eigs - plot computed eigenvalues
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
|
.seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
|
@*/
|
@*/
|
PetscErrorCode EPSSolve(EPS eps)
|
PetscErrorCode EPSSolve(EPS eps)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
PetscReal re,im;
|
PetscReal re,im;
|
PetscTruth flg;
|
PetscTruth flg;
|
PetscViewer viewer;
|
PetscViewer viewer;
|
PetscDraw draw;
|
PetscDraw draw;
|
PetscDrawSP drawsp;
|
PetscDrawSP drawsp;
|
STMatMode matmode;
|
STMatMode matmode;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
|
|
ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_view_binary",&flg);CHKERRQ(ierr);
|
ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_view_binary",&flg);CHKERRQ(ierr);
|
if (flg) {
|
if (flg) {
|
Mat A,B;
|
Mat A,B;
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = MatView(A,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));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);
|
if (B) ierr = MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)eps)->comm));CHKERRQ(ierr);
|
}
|
}
|
|
|
/* reset the convergence flag from the previous solves */
|
/* reset the convergence flag from the previous solves */
|
eps->reason = EPS_CONVERGED_ITERATING;
|
eps->reason = EPS_CONVERGED_ITERATING;
|
|
|
if (!eps->setupcalled){ ierr = EPSSetUp(eps);CHKERRQ(ierr); }
|
if (!eps->setupcalled){ ierr = EPSSetUp(eps);CHKERRQ(ierr); }
|
ierr = STResetOperationCounters(eps->OP);CHKERRQ(ierr);
|
ierr = STResetOperationCounters(eps->OP);CHKERRQ(ierr);
|
ierr = IPResetOperationCounters(eps->ip);CHKERRQ(ierr);
|
ierr = IPResetOperationCounters(eps->ip);CHKERRQ(ierr);
|
eps->evecsavailable = PETSC_FALSE;
|
eps->evecsavailable = PETSC_FALSE;
|
eps->nconv = 0;
|
eps->nconv = 0;
|
eps->its = 0;
|
eps->its = 0;
|
for (i=0;i<eps->ncv;i++) eps->eigr[i]=eps->eigi[i]=eps->errest[i]=0.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);
|
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);
|
ierr = PetscLogEventBegin(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);
|
|
|
switch (eps->solverclass) {
|
switch (eps->solverclass) {
|
case EPS_ONE_SIDE:
|
case EPS_ONE_SIDE:
|
ierr = (*eps->ops->solve)(eps);CHKERRQ(ierr); break;
|
ierr = (*eps->ops->solve)(eps);CHKERRQ(ierr); break;
|
case EPS_TWO_SIDE:
|
case EPS_TWO_SIDE:
|
if (eps->ops->solvets) {
|
if (eps->ops->solvets) {
|
ierr = (*eps->ops->solvets)(eps);CHKERRQ(ierr); break;
|
ierr = (*eps->ops->solvets)(eps);CHKERRQ(ierr); break;
|
} else {
|
} else {
|
SETERRQ(1,"Two-sided version unavailable for this solver");
|
SETERRQ(1,"Two-sided version unavailable for this solver");
|
}
|
}
|
default:
|
default:
|
SETERRQ(1,"Wrong value of eps->solverclass");
|
SETERRQ(1,"Wrong value of eps->solverclass");
|
}
|
}
|
|
|
ierr = STGetMatMode(eps->OP,&matmode);CHKERRQ(ierr);
|
ierr = STGetMatMode(eps->OP,&matmode);CHKERRQ(ierr);
|
if (matmode == STMATMODE_INPLACE && eps->ispositive) {
|
if (matmode == STMATMODE_INPLACE && eps->ispositive) {
|
/* Purge eigenvectors before reverting operator */
|
/* Purge eigenvectors before reverting operator */
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
}
|
}
|
ierr = STPostSolve(eps->OP);CHKERRQ(ierr);
|
ierr = STPostSolve(eps->OP);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Solve,eps,eps->V[0],eps->V[0],0);CHKERRQ(ierr);
|
|
|
if (!eps->reason) {
|
if (!eps->reason) {
|
SETERRQ(1,"Internal error, solver returned without setting converged reason");
|
SETERRQ(1,"Internal error, solver returned without setting converged reason");
|
}
|
}
|
|
|
/* Map eigenvalues back to the original problem, necessary in some
|
/* Map eigenvalues back to the original problem, necessary in some
|
* spectral transformations */
|
* spectral transformations */
|
ierr = (*eps->ops->backtransform)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->backtransform)(eps);CHKERRQ(ierr);
|
|
|
/* Adjust left eigenvectors in generalized problems: y = B^T y */
|
/* Adjust left eigenvectors in generalized problems: y = B^T y */
|
if (eps->isgeneralized && eps->solverclass == EPS_TWO_SIDE) {
|
if (eps->isgeneralized && eps->solverclass == EPS_TWO_SIDE) {
|
Mat B;
|
Mat B;
|
KSP ksp;
|
KSP ksp;
|
Vec w;
|
Vec w;
|
ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);
|
ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);
|
ierr = KSPCreate(((PetscObject)eps)->comm,&ksp);CHKERRQ(ierr);
|
ierr = KSPCreate(((PetscObject)eps)->comm,&ksp);CHKERRQ(ierr);
|
ierr = KSPSetOperators(ksp,B,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
|
ierr = KSPSetOperators(ksp,B,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
|
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
|
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
|
ierr = MatGetVecs(B,PETSC_NULL,&w);CHKERRQ(ierr);
|
ierr = MatGetVecs(B,PETSC_NULL,&w);CHKERRQ(ierr);
|
for (i=0;i<eps->nconv;i++) {
|
for (i=0;i<eps->nconv;i++) {
|
ierr = VecCopy(eps->W[i],w);CHKERRQ(ierr);
|
ierr = VecCopy(eps->W[i],w);CHKERRQ(ierr);
|
ierr = KSPSolveTranspose(ksp,w,eps->W[i]);CHKERRQ(ierr);
|
ierr = KSPSolveTranspose(ksp,w,eps->W[i]);CHKERRQ(ierr);
|
}
|
}
|
ierr = KSPDestroy(ksp);CHKERRQ(ierr);
|
ierr = KSPDestroy(ksp);CHKERRQ(ierr);
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
}
|
}
|
|
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
/* reorder conjugate eigenvalues (positive imaginary first) */
|
/* reorder conjugate eigenvalues (positive imaginary first) */
|
for (i=0; i<eps->nconv-1; i++) {
|
for (i=0; i<eps->nconv-1; i++) {
|
if (eps->eigi[i] != 0) {
|
if (eps->eigi[i] != 0) {
|
if (eps->eigi[i] < 0) {
|
if (eps->eigi[i] < 0) {
|
eps->eigi[i] = -eps->eigi[i];
|
eps->eigi[i] = -eps->eigi[i];
|
eps->eigi[i+1] = -eps->eigi[i+1];
|
eps->eigi[i+1] = -eps->eigi[i+1];
|
ierr = VecScale(eps->V[i+1],-1.0); CHKERRQ(ierr);
|
ierr = VecScale(eps->V[i+1],-1.0); CHKERRQ(ierr);
|
}
|
}
|
i++;
|
i++;
|
}
|
}
|
}
|
}
|
#endif
|
#endif
|
|
|
/* sort eigenvalues according to eps->which parameter */
|
/* sort eigenvalues according to eps->which parameter */
|
ierr = PetscFree(eps->perm);CHKERRQ(ierr);
|
ierr = PetscFree(eps->perm);CHKERRQ(ierr);
|
if (eps->nconv > 0) {
|
if (eps->nconv > 0) {
|
ierr = PetscMalloc(sizeof(PetscInt)*eps->nconv, &eps->perm); CHKERRQ(ierr);
|
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 = 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);
|
ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_view",&flg);CHKERRQ(ierr);
|
if (flg && !PetscPreLoadingOn) { ierr = EPSView(eps,PETSC_VIEWER_STDOUT_WORLD);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);
|
ierr = PetscOptionsHasName(((PetscObject)eps)->prefix,"-eps_plot_eigs",&flg);CHKERRQ(ierr);
|
if (flg) {
|
if (flg) {
|
ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Computed Eigenvalues",
|
ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Computed Eigenvalues",
|
PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);CHKERRQ(ierr);
|
PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);CHKERRQ(ierr);
|
ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
|
ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
|
ierr = PetscDrawSPCreate(draw,1,&drawsp);CHKERRQ(ierr);
|
ierr = PetscDrawSPCreate(draw,1,&drawsp);CHKERRQ(ierr);
|
for( i=0; i<eps->nconv; i++ ) {
|
for( i=0; i<eps->nconv; i++ ) {
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
re = PetscRealPart(eps->eigr[i]);
|
re = PetscRealPart(eps->eigr[i]);
|
im = PetscImaginaryPart(eps->eigi[i]);
|
im = PetscImaginaryPart(eps->eigi[i]);
|
#else
|
#else
|
re = eps->eigr[i];
|
re = eps->eigr[i];
|
im = eps->eigi[i];
|
im = eps->eigi[i];
|
#endif
|
#endif
|
ierr = PetscDrawSPAddPoint(drawsp,&re,&im);CHKERRQ(ierr);
|
ierr = PetscDrawSPAddPoint(drawsp,&re,&im);CHKERRQ(ierr);
|
}
|
}
|
ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr);
|
ierr = PetscDrawSPDraw(drawsp);CHKERRQ(ierr);
|
ierr = PetscDrawSPDestroy(drawsp);CHKERRQ(ierr);
|
ierr = PetscDrawSPDestroy(drawsp);CHKERRQ(ierr);
|
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
|
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
|
}
|
}
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetIterationNumber"
|
#define __FUNCT__ "EPSGetIterationNumber"
|
/*@
|
/*@
|
EPSGetIterationNumber - Gets the current iteration number. If the
|
EPSGetIterationNumber - Gets the current iteration number. If the
|
call to EPSSolve() is complete, then it returns the number of iterations
|
call to EPSSolve() is complete, then it returns the number of iterations
|
carried out by the solution method.
|
carried out by the solution method.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
|
|
Output Parameter:
|
Output Parameter:
|
. its - number of iterations
|
. its - number of iterations
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
Note:
|
Note:
|
During the i-th iteration this call returns i-1. If EPSSolve() is
|
During the i-th iteration this call returns i-1. If EPSSolve() is
|
complete, then parameter "its" contains either the iteration number at
|
complete, then parameter "its" contains either the iteration number at
|
which convergence was successfully reached, or failure was detected.
|
which convergence was successfully reached, or failure was detected.
|
Call EPSGetConvergedReason() to determine if the solver converged or
|
Call EPSGetConvergedReason() to determine if the solver converged or
|
failed and why.
|
failed and why.
|
|
|
.seealso: EPSGetConvergedReason(), EPSSetTolerances()
|
.seealso: EPSGetConvergedReason(), EPSSetTolerances()
|
@*/
|
@*/
|
PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
|
PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidIntPointer(its,2);
|
PetscValidIntPointer(its,2);
|
*its = eps->its;
|
*its = eps->its;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetOperationCounters"
|
#define __FUNCT__ "EPSGetOperationCounters"
|
/*@
|
/*@
|
EPSGetOperationCounters - Gets the total number of operator applications,
|
EPSGetOperationCounters - Gets the total number of operator applications,
|
inner product operations and linear iterations used by the ST object
|
inner product operations and linear iterations used by the ST object
|
during the last EPSSolve() call.
|
during the last EPSSolve() call.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - EPS context
|
. eps - EPS context
|
|
|
Output Parameter:
|
Output Parameter:
|
+ ops - number of operator applications
|
+ ops - number of operator applications
|
. dots - number of inner product operations
|
. dots - number of inner product operations
|
- lits - number of linear iterations
|
- lits - number of linear iterations
|
|
|
Notes:
|
Notes:
|
When the eigensolver algorithm invokes STApply() then a linear system
|
When the eigensolver algorithm invokes STApply() then a linear system
|
must be solved (except in the case of standard eigenproblems and shift
|
must be solved (except in the case of standard eigenproblems and shift
|
transformation). The number of iterations required in this solve is
|
transformation). The number of iterations required in this solve is
|
accumulated into a counter whose value is returned by this function.
|
accumulated into a counter whose value is returned by this function.
|
|
|
These counters are reset to zero at each successive call to EPSSolve().
|
These counters are reset to zero at each successive call to EPSSolve().
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
@*/
|
@*/
|
PetscErrorCode EPSGetOperationCounters(EPS eps,PetscInt* ops,PetscInt* dots,PetscInt* lits)
|
PetscErrorCode EPSGetOperationCounters(EPS eps,PetscInt* ops,PetscInt* dots,PetscInt* lits)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
ierr = STGetOperationCounters(eps->OP,ops,lits);CHKERRQ(ierr);
|
ierr = STGetOperationCounters(eps->OP,ops,lits);CHKERRQ(ierr);
|
if (dots) {
|
if (dots) {
|
ierr = IPGetOperationCounters(eps->ip,dots);CHKERRQ(ierr);
|
ierr = IPGetOperationCounters(eps->ip,dots);CHKERRQ(ierr);
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetConverged"
|
#define __FUNCT__ "EPSGetConverged"
|
/*@
|
/*@
|
EPSGetConverged - Gets the number of converged eigenpairs.
|
EPSGetConverged - Gets the number of converged eigenpairs.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
|
|
Output Parameter:
|
Output Parameter:
|
. nconv - number of converged eigenpairs
|
. nconv - number of converged eigenpairs
|
|
|
Note:
|
Note:
|
This function should be called after EPSSolve() has finished.
|
This function should be called after EPSSolve() has finished.
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSetDimensions()
|
.seealso: EPSSetDimensions()
|
@*/
|
@*/
|
PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
|
PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidIntPointer(nconv,2);
|
PetscValidIntPointer(nconv,2);
|
*nconv = eps->nconv;
|
*nconv = eps->nconv;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetConvergedReason"
|
#define __FUNCT__ "EPSGetConvergedReason"
|
/*@C
|
/*@C
|
EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
|
EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
|
stopped.
|
stopped.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
|
|
Output Parameter:
|
Output Parameter:
|
. reason - negative value indicates diverged, positive value converged
|
. reason - negative value indicates diverged, positive value converged
|
(see EPSConvergedReason)
|
(see EPSConvergedReason)
|
|
|
Possible values for reason:
|
Possible values for reason:
|
+ EPS_CONVERGED_TOL - converged up to tolerance
|
+ EPS_CONVERGED_TOL - converged up to tolerance
|
. EPS_DIVERGED_ITS - required more than its to reach convergence
|
. EPS_DIVERGED_ITS - required more than its to reach convergence
|
. EPS_DIVERGED_BREAKDOWN - generic breakdown in method
|
. EPS_DIVERGED_BREAKDOWN - generic breakdown in method
|
- EPS_DIVERGED_NONSYMMETRIC - The operator is nonsymmetric
|
- EPS_DIVERGED_NONSYMMETRIC - The operator is nonsymmetric
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
Notes: Can only be called after the call to EPSSolve() is complete.
|
Notes: Can only be called after the call to EPSSolve() is complete.
|
|
|
.seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
|
.seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
|
@*/
|
@*/
|
PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
|
PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidIntPointer(reason,2);
|
PetscValidIntPointer(reason,2);
|
*reason = eps->reason;
|
*reason = eps->reason;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetInvariantSubspace"
|
#define __FUNCT__ "EPSGetInvariantSubspace"
|
/*@
|
/*@
|
EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
|
EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
|
subspace.
|
subspace.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
|
|
Output Parameter:
|
Output Parameter:
|
. v - an array of vectors
|
. v - an array of vectors
|
|
|
Notes:
|
Notes:
|
This function should be called after EPSSolve() has finished.
|
This function should be called after EPSSolve() has finished.
|
|
|
The user should provide in v an array of nconv vectors, where nconv is
|
The user should provide in v an array of nconv vectors, where nconv is
|
the value returned by EPSGetConverged().
|
the value returned by EPSGetConverged().
|
|
|
The first k vectors returned in v span an invariant subspace associated
|
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
|
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-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
|
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
|
in X for all x in X (a similar definition applies for generalized
|
eigenproblems).
|
eigenproblems).
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetLeftInvariantSubspace()
|
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetLeftInvariantSubspace()
|
@*/
|
@*/
|
PetscErrorCode EPSGetInvariantSubspace(EPS eps, Vec *v)
|
PetscErrorCode EPSGetInvariantSubspace(EPS eps, Vec *v)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidPointer(v,2);
|
PetscValidPointer(v,2);
|
PetscValidHeaderSpecific(*v,VEC_COOKIE,2);
|
PetscValidHeaderSpecific(*v,VEC_COOKIE,2);
|
if (!eps->V) {
|
if (!eps->V) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (!eps->ishermitian && eps->evecsavailable) {
|
if (!eps->ishermitian && eps->evecsavailable) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetRightVector,EPSComputeRelativeError or EPSComputeResidualNorm");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetRightVector,EPSComputeRelativeError or EPSComputeResidualNorm");
|
}
|
}
|
for (i=0;i<eps->nconv;i++) {
|
for (i=0;i<eps->nconv;i++) {
|
ierr = VecCopy(eps->V[i],v[i]);CHKERRQ(ierr);
|
ierr = VecCopy(eps->V[i],v[i]);CHKERRQ(ierr);
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetLeftInvariantSubspace"
|
#define __FUNCT__ "EPSGetLeftInvariantSubspace"
|
/*@
|
/*@
|
EPSGetLeftInvariantSubspace - Gets an orthonormal basis of the computed left
|
EPSGetLeftInvariantSubspace - Gets an orthonormal basis of the computed left
|
invariant subspace (only available in two-sided eigensolvers).
|
invariant subspace (only available in two-sided eigensolvers).
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
|
|
Output Parameter:
|
Output Parameter:
|
. v - an array of vectors
|
. v - an array of vectors
|
|
|
Notes:
|
Notes:
|
This function should be called after EPSSolve() has finished.
|
This function should be called after EPSSolve() has finished.
|
|
|
The user should provide in v an array of nconv vectors, where nconv is
|
The user should provide in v an array of nconv vectors, where nconv is
|
the value returned by EPSGetConverged().
|
the value returned by EPSGetConverged().
|
|
|
The first k vectors returned in v span a left invariant subspace associated
|
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
|
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-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
|
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
|
in Y for all y in Y (a similar definition applies for generalized
|
eigenproblems).
|
eigenproblems).
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetInvariantSubspace
|
.seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve(), EPSGetInvariantSubspace
|
@*/
|
@*/
|
PetscErrorCode EPSGetLeftInvariantSubspace(EPS eps, Vec *v)
|
PetscErrorCode EPSGetLeftInvariantSubspace(EPS eps, Vec *v)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidPointer(v,2);
|
PetscValidPointer(v,2);
|
PetscValidHeaderSpecific(*v,VEC_COOKIE,2);
|
PetscValidHeaderSpecific(*v,VEC_COOKIE,2);
|
if (!eps->W) {
|
if (!eps->W) {
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
} else {
|
} else {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
}
|
}
|
for (i=0;i<eps->nconv;i++) {
|
for (i=0;i<eps->nconv;i++) {
|
ierr = VecCopy(eps->W[i],v[i]);CHKERRQ(ierr);
|
ierr = VecCopy(eps->W[i],v[i]);CHKERRQ(ierr);
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetEigenpair"
|
#define __FUNCT__ "EPSGetEigenpair"
|
/*@
|
/*@
|
EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
|
EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
|
EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
|
EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of the solution
|
- i - index of the solution
|
|
|
Output Parameters:
|
Output Parameters:
|
+ eigr - real part of eigenvalue
|
+ eigr - real part of eigenvalue
|
. eigi - imaginary part of eigenvalue
|
. eigi - imaginary part of eigenvalue
|
. Vr - real part of eigenvector
|
. Vr - real part of eigenvector
|
- Vi - imaginary part of eigenvector
|
- Vi - imaginary part of eigenvector
|
|
|
Notes:
|
Notes:
|
If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
|
If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
|
configured with complex scalars the eigenvalue is stored
|
configured with complex scalars the eigenvalue is stored
|
directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
|
directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
|
set to zero).
|
set to zero).
|
|
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSGetValue(), EPSGetRightVector(), EPSGetLeftVector(), EPSSolve(),
|
.seealso: EPSGetValue(), EPSGetRightVector(), EPSGetLeftVector(), EPSSolve(),
|
EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
|
EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
|
@*/
|
@*/
|
PetscErrorCode EPSGetEigenpair(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi, Vec Vr, Vec Vi)
|
PetscErrorCode EPSGetEigenpair(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi, Vec Vr, Vec Vi)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->eigr || !eps->eigi || !eps->V) {
|
if (!eps->eigr || !eps->eigi || !eps->V) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
ierr = EPSGetValue(eps,i,eigr,eigi);CHKERRQ(ierr);
|
ierr = EPSGetValue(eps,i,eigr,eigi);CHKERRQ(ierr);
|
ierr = EPSGetRightVector(eps,i,Vr,Vi);CHKERRQ(ierr);
|
ierr = EPSGetRightVector(eps,i,Vr,Vi);CHKERRQ(ierr);
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetValue"
|
#define __FUNCT__ "EPSGetValue"
|
/*@
|
/*@
|
EPSGetValue - Gets the i-th eigenvalue as computed by EPSSolve().
|
EPSGetValue - Gets the i-th eigenvalue as computed by EPSSolve().
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of the solution
|
- i - index of the solution
|
|
|
Output Parameters:
|
Output Parameters:
|
+ eigr - real part of eigenvalue
|
+ eigr - real part of eigenvalue
|
- eigi - imaginary part of eigenvalue
|
- eigi - imaginary part of eigenvalue
|
|
|
Notes:
|
Notes:
|
If the eigenvalue is real, then eigi is set to zero. If PETSc is
|
If the eigenvalue is real, then eigi is set to zero. If PETSc is
|
configured with complex scalars the eigenvalue is stored
|
configured with complex scalars the eigenvalue is stored
|
directly in eigr (eigi is set to zero).
|
directly in eigr (eigi is set to zero).
|
|
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
EPSGetEigenpair()
|
EPSGetEigenpair()
|
@*/
|
@*/
|
PetscErrorCode EPSGetValue(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi)
|
PetscErrorCode EPSGetValue(EPS eps, PetscInt i, PetscScalar *eigr, PetscScalar *eigi)
|
{
|
{
|
PetscInt k;
|
PetscInt k;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->eigr || !eps->eigi) {
|
if (!eps->eigr || !eps->eigi) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
|
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#ifdef PETSC_USE_COMPLEX
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigi) *eigi = 0;
|
if (eigi) *eigi = 0;
|
#else
|
#else
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigi) *eigi = eps->eigi[k];
|
if (eigi) *eigi = eps->eigi[k];
|
#endif
|
#endif
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetRightVector"
|
#define __FUNCT__ "EPSGetRightVector"
|
/*@
|
/*@
|
EPSGetRightVector - Gets the i-th right eigenvector as computed by EPSSolve().
|
EPSGetRightVector - Gets the i-th right eigenvector as computed by EPSSolve().
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of the solution
|
- i - index of the solution
|
|
|
Output Parameters:
|
Output Parameters:
|
+ Vr - real part of eigenvector
|
+ Vr - real part of eigenvector
|
- Vi - imaginary part of eigenvector
|
- Vi - imaginary part of eigenvector
|
|
|
Notes:
|
Notes:
|
If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
|
If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
|
configured with complex scalars the eigenvector is stored
|
configured with complex scalars the eigenvector is stored
|
directly in Vr (Vi is set to zero).
|
directly in Vr (Vi is set to zero).
|
|
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
EPSGetEigenpair(), EPSGetLeftVector()
|
EPSGetEigenpair(), EPSGetLeftVector()
|
@*/
|
@*/
|
PetscErrorCode EPSGetRightVector(EPS eps, PetscInt i, Vec Vr, Vec Vi)
|
PetscErrorCode EPSGetRightVector(EPS eps, PetscInt i, Vec Vr, Vec Vi)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt k;
|
PetscInt k;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->V) {
|
if (!eps->V) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
if (!eps->evecsavailable && (Vr || Vi) ) {
|
if (!eps->evecsavailable && (Vr || Vi) ) {
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
}
|
}
|
|
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#ifdef PETSC_USE_COMPLEX
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
#else
|
#else
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecCopy(eps->V[k+1], Vi); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecCopy(eps->V[k+1], Vi); CHKERRQ(ierr); }
|
} else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
|
} else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
|
if (Vr) { ierr = VecCopy(eps->V[k-1], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k-1], Vr); CHKERRQ(ierr); }
|
if (Vi) {
|
if (Vi) {
|
ierr = VecCopy(eps->V[k], Vi); CHKERRQ(ierr);
|
ierr = VecCopy(eps->V[k], Vi); CHKERRQ(ierr);
|
ierr = VecScale(Vi,-1.0); CHKERRQ(ierr);
|
ierr = VecScale(Vi,-1.0); CHKERRQ(ierr);
|
}
|
}
|
} else { /* real eigenvalue */
|
} else { /* real eigenvalue */
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
}
|
}
|
#endif
|
#endif
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetLeftVector"
|
#define __FUNCT__ "EPSGetLeftVector"
|
/*@
|
/*@
|
EPSGetLeftVector - Gets the i-th left eigenvector as computed by EPSSolve()
|
EPSGetLeftVector - Gets the i-th left eigenvector as computed by EPSSolve()
|
(only available in two-sided eigensolvers).
|
(only available in two-sided eigensolvers).
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of the solution
|
- i - index of the solution
|
|
|
Output Parameters:
|
Output Parameters:
|
+ Wr - real part of eigenvector
|
+ Wr - real part of eigenvector
|
- Wi - imaginary part of eigenvector
|
- Wi - imaginary part of eigenvector
|
|
|
Notes:
|
Notes:
|
If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
|
If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
|
configured with complex scalars the eigenvector is stored
|
configured with complex scalars the eigenvector is stored
|
directly in Wr (Wi is set to zero).
|
directly in Wr (Wi is set to zero).
|
|
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(),
|
EPSGetEigenpair(), EPSGetLeftVector()
|
EPSGetEigenpair(), EPSGetLeftVector()
|
@*/
|
@*/
|
PetscErrorCode EPSGetLeftVector(EPS eps, PetscInt i, Vec Wr, Vec Wi)
|
PetscErrorCode EPSGetLeftVector(EPS eps, PetscInt i, Vec Wr, Vec Wi)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt k;
|
PetscInt k;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->W) {
|
if (!eps->W) {
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
} else {
|
} else {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
if (!eps->evecsavailable && (Wr || Wi) ) {
|
if (!eps->evecsavailable && (Wr || Wi) ) {
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
}
|
}
|
|
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#ifdef PETSC_USE_COMPLEX
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
#else
|
#else
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecCopy(eps->W[k+1], Wi); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecCopy(eps->W[k+1], Wi); CHKERRQ(ierr); }
|
} else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
|
} else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
|
if (Wr) { ierr = VecCopy(eps->W[k-1], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k-1], Wr); CHKERRQ(ierr); }
|
if (Wi) {
|
if (Wi) {
|
ierr = VecCopy(eps->W[k], Wi); CHKERRQ(ierr);
|
ierr = VecCopy(eps->W[k], Wi); CHKERRQ(ierr);
|
ierr = VecScale(Wi,-1.0); CHKERRQ(ierr);
|
ierr = VecScale(Wi,-1.0); CHKERRQ(ierr);
|
}
|
}
|
} else { /* real eigenvalue */
|
} else { /* real eigenvalue */
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
}
|
}
|
#endif
|
#endif
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetErrorEstimate"
|
#define __FUNCT__ "EPSGetErrorEstimate"
|
/*@
|
/*@
|
EPSGetErrorEstimate - Returns the error estimate associated to the i-th
|
EPSGetErrorEstimate - Returns the error estimate associated to the i-th
|
computed eigenpair.
|
computed eigenpair.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of eigenpair
|
- i - index of eigenpair
|
|
|
Output Parameter:
|
Output Parameter:
|
. errest - the error estimate
|
. errest - the error estimate
|
|
|
Notes:
|
Notes:
|
This is the error estimate used internally by the eigensolver. The actual
|
This is the error estimate used internally by the eigensolver. The actual
|
error bound can be computed with EPSComputeRelativeError(). See also the user's
|
error bound can be computed with EPSComputeRelativeError(). See also the user's
|
manual for details.
|
manual for details.
|
|
|
Level: advanced
|
Level: advanced
|
|
|
.seealso: EPSComputeRelativeError()
|
.seealso: EPSComputeRelativeError()
|
@*/
|
@*/
|
PetscErrorCode EPSGetErrorEstimate(EPS eps, PetscInt i, PetscReal *errest)
|
PetscErrorCode EPSGetErrorEstimate(EPS eps, PetscInt i, PetscReal *errest)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->eigr || !eps->eigi) {
|
if (!eps->eigr || !eps->eigi) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
if (eps->perm) i = eps->perm[i];
|
if (eps->perm) i = eps->perm[i];
|
if (errest) *errest = eps->errest[i];
|
if (errest) *errest = eps->errest[i];
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetErrorEstimateLeft"
|
#define __FUNCT__ "EPSGetErrorEstimateLeft"
|
/*@
|
/*@
|
EPSGetErrorEstimateLeft - Returns the left error estimate associated to the i-th
|
EPSGetErrorEstimateLeft - Returns the left error estimate associated to the i-th
|
computed eigenpair (only available in two-sided eigensolvers).
|
computed eigenpair (only available in two-sided eigensolvers).
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
+ eps - eigensolver context
|
+ eps - eigensolver context
|
- i - index of eigenpair
|
- i - index of eigenpair
|
|
|
Output Parameter:
|
Output Parameter:
|
. errest - the left error estimate
|
. errest - the left error estimate
|
|
|
Notes:
|
Notes:
|
This is the error estimate used internally by the eigensolver. The actual
|
This is the error estimate used internally by the eigensolver. The actual
|
error bound can be computed with EPSComputeRelativeErrorLeft(). See also the user's
|
error bound can be computed with EPSComputeRelativeErrorLeft(). See also the user's
|
manual for details.
|
manual for details.
|
|
|
Level: advanced
|
Level: advanced
|
|
|
.seealso: EPSComputeRelativeErrorLeft()
|
.seealso: EPSComputeRelativeErrorLeft()
|
@*/
|
@*/
|
PetscErrorCode EPSGetErrorEstimateLeft(EPS eps, PetscInt i, PetscReal *errest)
|
PetscErrorCode EPSGetErrorEstimateLeft(EPS eps, PetscInt i, PetscReal *errest)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
if (!eps->eigr || !eps->eigi) {
|
if (!eps->eigr || !eps->eigi) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "EPSSolve must be called first");
|
}
|
}
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
if (eps->solverclass!=EPS_TWO_SIDE) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Only available for two-sided solvers");
|
}
|
}
|
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
if (eps->perm) i = eps->perm[i];
|
if (eps->perm) i = eps->perm[i];
|
if (errest) *errest = eps->errest_left[i];
|
if (errest) *errest = eps->errest_left[i];
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSComputeResidualNorm"
|
#define __FUNCT__ "EPSComputeResidualNorm"
|
/*@
|
/*@
|
EPSComputeResidualNorm - Computes the norm of the residual vector associated with
|
EPSComputeResidualNorm - Computes the norm of the residual vector associated with
|
the i-th computed eigenpair.
|
the i-th computed eigenpair.
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
. i - the solution index
|
. i - the solution index
|
|
|
Output Parameter:
|
Output Parameter:
|
. norm - the residual norm, computed as ||Ax-kBx||_2 where k is the
|
. norm - the residual norm, computed as ||Ax-kBx||_2 where k is the
|
eigenvalue and x is the eigenvector.
|
eigenvalue and x is the eigenvector.
|
If k=0 then the residual norm is computed as ||Ax||_2.
|
If k=0 then the residual norm is computed as ||Ax||_2.
|
|
|
Notes:
|
Notes:
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
|
@*/
|
@*/
|
PetscErrorCode EPSComputeResidualNorm(EPS eps, PetscInt i, PetscReal *norm)
|
PetscErrorCode EPSComputeResidualNorm(EPS eps, PetscInt i, PetscReal *norm)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec u, v, w, xr, xi;
|
Vec u, v, w, xr, xi;
|
Mat A, B;
|
Mat A, B;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
PetscReal ni, nr;
|
PetscReal ni, nr;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&u); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&u); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&v); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&v); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&w); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&w); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = MatMult( A, xr, u ); CHKERRQ(ierr); /* u=A*x */
|
ierr = MatMult( A, xr, u ); CHKERRQ(ierr); /* u=A*x */
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (eps->isgeneralized) { ierr = MatMult( B, xr, w ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMult( B, xr, w ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B*x */
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B*x */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*x-k*B*x */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*x-k*B*x */
|
}
|
}
|
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
} else {
|
} else {
|
ierr = MatMult( A, xr, u ); CHKERRQ(ierr); /* u=A*xr */
|
ierr = MatMult( A, xr, u ); CHKERRQ(ierr); /* u=A*xr */
|
if (eps->isgeneralized) { ierr = MatMult( B, xr, v ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMult( B, xr, v ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B*xr */
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr */
|
if (eps->isgeneralized) { ierr = MatMult( B, xi, w ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMult( B, xi, w ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B*xi */
|
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B*xi */
|
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr+ki*B*xi */
|
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr+ki*B*xi */
|
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
|
ierr = MatMult( A, xi, u ); CHKERRQ(ierr); /* u=A*xi */
|
ierr = MatMult( A, xi, u ); CHKERRQ(ierr); /* u=A*xi */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi */
|
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi-ki*B*xr */
|
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi-ki*B*xr */
|
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
|
*norm = SlepcAbsEigenvalue( nr, ni );
|
*norm = SlepcAbsEigenvalue( nr, ni );
|
}
|
}
|
#endif
|
#endif
|
|
|
ierr = VecDestroy(w); CHKERRQ(ierr);
|
ierr = VecDestroy(w); CHKERRQ(ierr);
|
ierr = VecDestroy(v); CHKERRQ(ierr);
|
ierr = VecDestroy(v); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSComputeResidualNormLeft"
|
#define __FUNCT__ "EPSComputeResidualNormLeft"
|
/*@
|
/*@
|
EPSComputeResidualNormLeft - Computes the norm of the residual vector associated with
|
EPSComputeResidualNormLeft - Computes the norm of the residual vector associated with
|
the i-th computed left eigenvector (only available in two-sided eigensolvers).
|
the i-th computed left eigenvector (only available in two-sided eigensolvers).
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
. i - the solution index
|
. i - the solution index
|
|
|
Output Parameter:
|
Output Parameter:
|
. norm - the residual norm, computed as ||y'A-ky'B||_2 where k is the
|
. norm - the residual norm, computed as ||y'A-ky'B||_2 where k is the
|
eigenvalue and y is the left eigenvector.
|
eigenvalue and y is the left eigenvector.
|
If k=0 then the residual norm is computed as ||y'A||_2.
|
If k=0 then the residual norm is computed as ||y'A||_2.
|
|
|
Notes:
|
Notes:
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
|
Eigenpairs are indexed according to the ordering criterion established
|
Eigenpairs are indexed according to the ordering criterion established
|
with EPSSetWhichEigenpairs().
|
with EPSSetWhichEigenpairs().
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
|
.seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs()
|
@*/
|
@*/
|
PetscErrorCode EPSComputeResidualNormLeft(EPS eps, PetscInt i, PetscReal *norm)
|
PetscErrorCode EPSComputeResidualNormLeft(EPS eps, PetscInt i, PetscReal *norm)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec u, v, w, xr, xi;
|
Vec u, v, w, xr, xi;
|
Mat A, B;
|
Mat A, B;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
PetscReal ni, nr;
|
PetscReal ni, nr;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&u); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&u); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&v); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&v); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&w); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&w); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xi); CHKERRQ(ierr);
|
ierr = EPSGetValue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetValue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetLeftVector(eps,i,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetLeftVector(eps,i,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*x */
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*x */
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, w ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, w ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B'*x */
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B'*x */
|
ierr = VecAXPY( u, -kr, w); CHKERRQ(ierr); /* u=A'*x-k*B'*x */
|
ierr = VecAXPY( u, -kr, w); CHKERRQ(ierr); /* u=A'*x-k*B'*x */
|
}
|
}
|
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
} else {
|
} else {
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*xr */
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*xr */
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, v ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, v ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B'*xr */
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B'*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xi, w ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xi, w ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B'*xi */
|
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B'*xi */
|
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr+ki*B'*xi */
|
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr+ki*B'*xi */
|
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
|
ierr = MatMultTranspose( A, xi, u ); CHKERRQ(ierr); /* u=A'*xi */
|
ierr = MatMultTranspose( A, xi, u ); CHKERRQ(ierr); /* u=A'*xi */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi */
|
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi */
|
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi-ki*B'*xr */
|
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi-ki*B'*xr */
|
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
|
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
|
*norm = SlepcAbsEigenvalue( nr, ni );
|
*norm = SlepcAbsEigenvalue( nr, ni );
|
}
|
}
|
#endif
|
#endif
|
|
|
ierr = VecDestroy(w); CHKERRQ(ierr);
|
ierr = VecDestroy(w); CHKERRQ(ierr);
|
ierr = VecDestroy(v); CHKERRQ(ierr);
|
ierr = VecDestroy(v); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSComputeRelativeError"
|
#define __FUNCT__ "EPSComputeRelativeError"
|
/*@
|
/*@
|
EPSComputeRelativeError - Computes the relative error bound associated
|
EPSComputeRelativeError - Computes the relative error bound associated
|
with the i-th computed eigenpair.
|
with the i-th computed eigenpair.
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
. i - the solution index
|
. i - the solution index
|
|
|
Output Parameter:
|
Output Parameter:
|
. error - the relative error bound, computed as ||Ax-kBx||_2/||kx||_2 where
|
. error - the relative error bound, computed as ||Ax-kBx||_2/||kx||_2 where
|
k is the eigenvalue and x is the eigenvector.
|
k is the eigenvalue and x is the eigenvector.
|
If k=0 the relative error is computed as ||Ax||_2/||x||_2.
|
If k=0 the relative error is computed as ||Ax||_2/||x||_2.
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSComputeResidualNorm(), EPSGetErrorEstimate()
|
.seealso: EPSSolve(), EPSComputeResidualNorm(), EPSGetErrorEstimate()
|
@*/
|
@*/
|
PetscErrorCode EPSComputeRelativeError(EPS eps, PetscInt i, PetscReal *error)
|
PetscErrorCode EPSComputeRelativeError(EPS eps, PetscInt i, PetscReal *error)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec xr, xi;
|
Vec xr, xi;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
PetscReal norm, er;
|
PetscReal norm, er;
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
Vec u;
|
Vec u;
|
PetscReal ei;
|
PetscReal ei;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
ierr = EPSComputeResidualNorm(eps,i,&norm); CHKERRQ(ierr);
|
ierr = EPSComputeResidualNorm(eps,i,&norm); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
if (PetscAbsScalar(kr) > norm) {
|
if (PetscAbsScalar(kr) > norm) {
|
*error = norm / (PetscAbsScalar(kr) * er);
|
*error = norm / (PetscAbsScalar(kr) * er);
|
} else {
|
} else {
|
*error = norm / er;
|
*error = norm / er;
|
}
|
}
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
} else {
|
} else {
|
if (SlepcAbsEigenvalue(kr,ki) > norm) {
|
if (SlepcAbsEigenvalue(kr,ki) > norm) {
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
|
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
|
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
|
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
} else {
|
} else {
|
ierr = VecDot(xr, xr, &er); CHKERRQ(ierr);
|
ierr = VecDot(xr, xr, &er); CHKERRQ(ierr);
|
ierr = VecDot(xi, xi, &ei); CHKERRQ(ierr);
|
ierr = VecDot(xi, xi, &ei); CHKERRQ(ierr);
|
}
|
}
|
*error = norm / SlepcAbsEigenvalue(er, ei);
|
*error = norm / SlepcAbsEigenvalue(er, ei);
|
}
|
}
|
#endif
|
#endif
|
|
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSComputeRelativeErrorLeft"
|
#define __FUNCT__ "EPSComputeRelativeErrorLeft"
|
/*@
|
/*@
|
EPSComputeRelativeErrorLeft - Computes the relative error bound associated
|
EPSComputeRelativeErrorLeft - Computes the relative error bound associated
|
with the i-th computed eigenvalue and left eigenvector (only available in
|
with the i-th computed eigenvalue and left eigenvector (only available in
|
two-sided eigensolvers).
|
two-sided eigensolvers).
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. eps - the eigensolver context
|
. eps - the eigensolver context
|
. i - the solution index
|
. i - the solution index
|
|
|
Output Parameter:
|
Output Parameter:
|
. error - the relative error bound, computed as ||y'A-ky'B||_2/||ky||_2 where
|
. error - the relative error bound, computed as ||y'A-ky'B||_2/||ky||_2 where
|
k is the eigenvalue and y is the left eigenvector.
|
k is the eigenvalue and y is the left eigenvector.
|
If k=0 the relative error is computed as ||y'A||_2/||y||_2.
|
If k=0 the relative error is computed as ||y'A||_2/||y||_2.
|
|
|
Level: beginner
|
Level: beginner
|
|
|
.seealso: EPSSolve(), EPSComputeResidualNormLeft(), EPSGetErrorEstimateLeft()
|
.seealso: EPSSolve(), EPSComputeResidualNormLeft(), EPSGetErrorEstimateLeft()
|
@*/
|
@*/
|
PetscErrorCode EPSComputeRelativeErrorLeft(EPS eps, PetscInt i, PetscReal *error)
|
PetscErrorCode EPSComputeRelativeErrorLeft(EPS eps, PetscInt i, PetscReal *error)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec xr, xi;
|
Vec xr, xi;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
PetscReal norm, er;
|
PetscReal norm, er;
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
Vec u;
|
Vec u;
|
PetscReal ei;
|
PetscReal ei;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
ierr = EPSComputeResidualNormLeft(eps,i,&norm); CHKERRQ(ierr);
|
ierr = EPSComputeResidualNormLeft(eps,i,&norm); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&xi); CHKERRQ(ierr);
|
ierr = EPSGetValue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetValue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetLeftVector(eps,i,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetLeftVector(eps,i,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
*error = norm / (PetscAbsScalar(kr) * er);
|
*error = norm / (PetscAbsScalar(kr) * er);
|
} else {
|
} else {
|
*error = norm / er;
|
*error = norm / er;
|
}
|
}
|
#ifndef PETSC_USE_COMPLEX
|
#ifndef PETSC_USE_COMPLEX
|
} else {
|
} else {
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
|
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
|
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
|
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
ierr = VecDestroy(u); CHKERRQ(ierr);
|
*error = norm / SlepcAbsEigenvalue(er, ei);
|
*error = norm / SlepcAbsEigenvalue(er, ei);
|
}
|
}
|
#endif
|
#endif
|
|
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xr); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
ierr = VecDestroy(xi); CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#define SWAP(a,b,t) {t=a;a=b;b=t;}
|
#define SWAP(a,b,t) {t=a;a=b;b=t;}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSSortEigenvalues"
|
#define __FUNCT__ "EPSSortEigenvalues"
|
/*@
|
/*@
|
EPSSortEigenvalues - Sorts a list of eigenvalues according to a certain
|
EPSSortEigenvalues - Sorts a list of eigenvalues according to a certain
|
criterion.
|
criterion.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - number of eigenvalue in the list
|
+ n - number of eigenvalue in the list
|
. eig - pointer to the array containing the eigenvalues
|
. eig - pointer to the array containing the eigenvalues
|
. eigi - imaginary part of the eigenvalues (only when using real numbers)
|
. eigi - imaginary part of the eigenvalues (only when using real numbers)
|
. which - sorting criterion
|
. which - sorting criterion
|
- nev - number of wanted eigenvalues
|
- nev - number of wanted eigenvalues
|
|
|
Output Parameter:
|
Output Parameter:
|
. permout - resulting permutation
|
. permout - resulting permutation
|
|
|
Notes:
|
Notes:
|
The result is a list of indices in the original eigenvalue array
|
The result is a list of indices in the original eigenvalue array
|
corresponding to the first nev eigenvalues sorted in the specified
|
corresponding to the first nev eigenvalues sorted in the specified
|
criterion
|
criterion
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSSortEigenvaluesReal(), EPSSetWhichEigenpairs()
|
.seealso: EPSSortEigenvaluesReal(), EPSSetWhichEigenpairs()
|
@*/
|
@*/
|
PetscErrorCode EPSSortEigenvalues(PetscInt n,PetscScalar *eig,PetscScalar *eigi,EPSWhich which,PetscInt nev,PetscInt *permout)
|
PetscErrorCode EPSSortEigenvalues(PetscInt n,PetscScalar *eig,PetscScalar *eigi,EPSWhich which,PetscInt nev,PetscInt *permout)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
PetscInt *perm;
|
PetscInt *perm;
|
PetscReal *values;
|
PetscReal *values;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&values);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&values);CHKERRQ(ierr);
|
for (i=0; i<n; i++) { perm[i] = i;}
|
for (i=0; i<n; i++) { perm[i] = i;}
|
|
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
for (i=0; i<n; i++) { values[i] = SlepcAbsEigenvalue(eig[i],eigi[i]); }
|
for (i=0; i<n; i++) { values[i] = SlepcAbsEigenvalue(eig[i],eigi[i]); }
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
for (i=0; i<n; i++) { values[i] = PetscRealPart(eig[i]); }
|
for (i=0; i<n; i++) { values[i] = PetscRealPart(eig[i]); }
|
break;
|
break;
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
for (i=0; i<n; i++) { values[i] = PetscImaginaryPart(eig[i]); }
|
for (i=0; i<n; i++) { values[i] = PetscImaginaryPart(eig[i]); }
|
#else
|
#else
|
for (i=0; i<n; i++) { values[i] = PetscAbsReal(eigi[i]); }
|
for (i=0; i<n; i++) { values[i] = PetscAbsReal(eigi[i]); }
|
#endif
|
#endif
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
|
|
ierr = PetscSortRealWithPermutation(n,values,perm);CHKERRQ(ierr);
|
ierr = PetscSortRealWithPermutation(n,values,perm);CHKERRQ(ierr);
|
|
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
for (i=0; i<nev; i++) { permout[i] = perm[n-1-i]; }
|
for (i=0; i<nev; i++) { permout[i] = perm[n-1-i]; }
|
break;
|
break;
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
for (i=0; i<nev; i++) { permout[i] = perm[i]; }
|
for (i=0; i<nev; i++) { permout[i] = perm[i]; }
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
|
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
for (i=0; i<nev-1; i++) {
|
for (i=0; i<nev-1; i++) {
|
if (eigi[permout[i]] != 0.0) {
|
if (eigi[permout[i]] != 0.0) {
|
if (eig[permout[i]] == eig[permout[i+1]] &&
|
if (eig[permout[i]] == eig[permout[i+1]] &&
|
eigi[permout[i]] == -eigi[permout[i+1]] &&
|
eigi[permout[i]] == -eigi[permout[i+1]] &&
|
eigi[permout[i]] < 0.0) {
|
eigi[permout[i]] < 0.0) {
|
PetscInt tmp;
|
PetscInt tmp;
|
SWAP(permout[i], permout[i+1], tmp);
|
SWAP(permout[i], permout[i+1], tmp);
|
}
|
}
|
i++;
|
i++;
|
}
|
}
|
}
|
}
|
#endif
|
#endif
|
|
|
ierr = PetscFree(values);CHKERRQ(ierr);
|
ierr = PetscFree(values);CHKERRQ(ierr);
|
ierr = PetscFree(perm);CHKERRQ(ierr);
|
ierr = PetscFree(perm);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSSortEigenvaluesReal"
|
#define __FUNCT__ "EPSSortEigenvaluesReal"
|
/*@
|
/*@
|
EPSSortEigenvaluesReal - Sorts a list of eigenvalues according to a certain
|
EPSSortEigenvaluesReal - Sorts a list of eigenvalues according to a certain
|
criterion (version for real eigenvalues only).
|
criterion (version for real eigenvalues only).
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - number of eigenvalue in the list
|
+ n - number of eigenvalue in the list
|
. eig - pointer to the array containing the eigenvalues (real)
|
. eig - pointer to the array containing the eigenvalues (real)
|
. which - sorting criterion
|
. which - sorting criterion
|
- nev - number of wanted eigenvalues
|
- nev - number of wanted eigenvalues
|
|
|
Output Parameter:
|
Output Parameter:
|
. permout - resulting permutation
|
. permout - resulting permutation
|
|
|
Workspace:
|
Workspace:
|
. work - workspace for storing n real values and n integer values
|
. work - workspace for storing n real values and n integer values
|
|
|
Notes:
|
Notes:
|
The result is a list of indices in the original eigenvalue array
|
The result is a list of indices in the original eigenvalue array
|
corresponding to the first nev eigenvalues sorted in the specified
|
corresponding to the first nev eigenvalues sorted in the specified
|
criterion
|
criterion
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSSortEigenvalues(), EPSSetWhichEigenpairs()
|
.seealso: EPSSortEigenvalues(), EPSSetWhichEigenpairs()
|
@*/
|
@*/
|
PetscErrorCode EPSSortEigenvaluesReal(PetscInt n,PetscReal *eig,EPSWhich which,PetscInt nev,PetscInt *permout,PetscReal *work)
|
PetscErrorCode EPSSortEigenvaluesReal(PetscInt n,PetscReal *eig,EPSWhich which,PetscInt nev,PetscInt *permout,PetscReal *work)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
PetscReal *values = work;
|
PetscReal *values = work;
|
PetscInt *perm = (PetscInt*)(work+n);
|
PetscInt *perm = (PetscInt*)(work+n);
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
for (i=0; i<n; i++) { perm[i] = i;}
|
for (i=0; i<n; i++) { perm[i] = i;}
|
|
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
for (i=0; i<n; i++) { values[i] = PetscAbsReal(eig[i]); }
|
for (i=0; i<n; i++) { values[i] = PetscAbsReal(eig[i]); }
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
for (i=0; i<n; i++) { values[i] = eig[i]; }
|
for (i=0; i<n; i++) { values[i] = eig[i]; }
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
|
|
ierr = PetscSortRealWithPermutation(n,values,perm);CHKERRQ(ierr);
|
ierr = PetscSortRealWithPermutation(n,values,perm);CHKERRQ(ierr);
|
|
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
for (i=0; i<nev; i++) { permout[i] = perm[n-1-i]; }
|
for (i=0; i<nev; i++) { permout[i] = perm[n-1-i]; }
|
break;
|
break;
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
for (i=0; i<nev; i++) { permout[i] = perm[i]; }
|
for (i=0; i<nev; i++) { permout[i] = perm[i]; }
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetStartVector"
|
#define __FUNCT__ "EPSGetStartVector"
|
/*@
|
/*@
|
EPSGetStartVector - Gets a vector to be used as the starting vector
|
EPSGetStartVector - Gets a vector to be used as the starting vector
|
in an Arnoldi or Lanczos reduction.
|
in an Arnoldi or Lanczos reduction.
|
|
|
Collective on EPS and Vec
|
Collective on EPS and Vec
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - the eigensolver context
|
+ eps - the eigensolver context
|
- i - index of the Arnoldi/Lanczos step
|
- i - index of the Arnoldi/Lanczos step
|
|
|
Output Parameters:
|
Output Parameters:
|
+ vec - the start vector
|
+ vec - the start vector
|
- breakdown - flag indicating that a breakdown has occurred
|
- breakdown - flag indicating that a breakdown has occurred
|
|
|
Notes:
|
Notes:
|
The start vector is computed from another vector: for the first step (i=0),
|
The start vector is computed from another vector: for the first step (i=0),
|
the initial vector is used (see EPSGetInitialVector()); otherwise a random
|
the initial vector is used (see EPSGetInitialVector()); otherwise a random
|
vector is created. Then this vector is forced to be in the range of OP (only
|
vector is created. Then this vector is forced to be in the range of OP (only
|
for generalized definite problems) and orthonormalized with respect to all
|
for generalized definite problems) and orthonormalized with respect to all
|
V-vectors up to i-1.
|
V-vectors up to i-1.
|
|
|
The flag breakdown is set to true if either i=0 and the vector belongs to the
|
The flag breakdown is set to true if either i=0 and the vector belongs to the
|
deflation space, or i>0 and the vector is linearly dependent with respect
|
deflation space, or i>0 and the vector is linearly dependent with respect
|
to the V-vectors.
|
to the V-vectors.
|
|
|
The caller must pass a vector already allocated with dimensions conforming
|
The caller must pass a vector already allocated with dimensions conforming
|
to the initial vector. This vector is overwritten.
|
to the initial vector. This vector is overwritten.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSGetInitialVector()
|
.seealso: EPSGetInitialVector()
|
|
|
@*/
|
@*/
|
PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,Vec vec,PetscTruth *breakdown)
|
PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,Vec vec,PetscTruth *breakdown)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal norm;
|
PetscReal norm;
|
PetscTruth lindep;
|
PetscTruth lindep;
|
Vec w;
|
Vec w;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(vec,VEC_COOKIE,3);
|
PetscValidHeaderSpecific(vec,VEC_COOKIE,3);
|
|
|
/* For the first step, use the initial vector, otherwise a random one */
|
/* For the first step, use the initial vector, otherwise a random one */
|
if (i==0) {
|
if (i==0) {
|
w = eps->vec_initial;
|
w = eps->vec_initial;
|
} else {
|
} else {
|
ierr = VecDuplicate(eps->vec_initial,&w);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial,&w);CHKERRQ(ierr);
|
ierr = SlepcVecSetRandom(w);CHKERRQ(ierr);
|
ierr = SlepcVecSetRandom(w);CHKERRQ(ierr);
|
}
|
}
|
|
|
/* Force the vector to be in the range of OP for definite generalized problems */
|
/* Force the vector to be in the range of OP for definite generalized problems */
|
if (eps->ispositive) {
|
if (eps->ispositive) {
|
ierr = STApply(eps->OP,w,vec);CHKERRQ(ierr);
|
ierr = STApply(eps->OP,w,vec);CHKERRQ(ierr);
|
} else {
|
} else {
|
ierr = VecCopy(w,vec);CHKERRQ(ierr);
|
ierr = VecCopy(w,vec);CHKERRQ(ierr);
|
}
|
}
|
|
|
/* Orthonormalize the vector with respect to previous vectors */
|
/* Orthonormalize the vector with respect to previous vectors */
|
ierr = IPOrthogonalize(eps->ip,i+eps->nds,PETSC_NULL,eps->DSV,vec,PETSC_NULL,&norm,&lindep,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(eps->ip,i+eps->nds,PETSC_NULL,eps->DSV,vec,PETSC_NULL,&norm,&lindep,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
if (breakdown) *breakdown = lindep;
|
if (breakdown) *breakdown = lindep;
|
else if (lindep || norm == 0.0) {
|
else if (lindep || norm == 0.0) {
|
if (i==0) { SETERRQ(1,"Initial vector is zero or belongs to the deflation space"); }
|
if (i==0) { SETERRQ(1,"Initial vector is zero or belongs to the deflation space"); }
|
else { SETERRQ(1,"Unable to generate more start vectors"); }
|
else { SETERRQ(1,"Unable to generate more start vectors"); }
|
}
|
}
|
|
|
ierr = VecScale(vec,1.0/norm);CHKERRQ(ierr);
|
ierr = VecScale(vec,1.0/norm);CHKERRQ(ierr);
|
|
|
if (i!=0) {
|
if (i!=0) {
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
}
|
}
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSGetLeftStartVector"
|
#define __FUNCT__ "EPSGetLeftStartVector"
|
/*@
|
/*@
|
EPSGetLeftStartVector - Gets a vector to be used as the starting vector
|
EPSGetLeftStartVector - Gets a vector to be used as the starting vector
|
in the left recurrence of a two-sided eigensolver.
|
in the left recurrence of a two-sided eigensolver.
|
|
|
Collective on EPS and Vec
|
Collective on EPS and Vec
|
|
|
Input Parameters:
|
Input Parameters:
|
+ eps - the eigensolver context
|
+ eps - the eigensolver context
|
- i - index of the Arnoldi/Lanczos step
|
- i - index of the Arnoldi/Lanczos step
|
|
|
Output Parameter:
|
Output Parameter:
|
. vec - the start vector
|
. vec - the start vector
|
|
|
Notes:
|
Notes:
|
The start vector is computed from another vector: for the first step (i=0),
|
The start vector is computed from another vector: for the first step (i=0),
|
the left initial vector is used (see EPSGetLeftInitialVector()); otherwise
|
the left initial vector is used (see EPSGetLeftInitialVector()); otherwise
|
a random vector is created. Then this vector is forced to be in the range
|
a random vector is created. Then this vector is forced to be in the range
|
of OP' and orthonormalized with respect to all W-vectors up to i-1.
|
of OP' and orthonormalized with respect to all W-vectors up to i-1.
|
|
|
The caller must pass a vector already allocated with dimensions conforming
|
The caller must pass a vector already allocated with dimensions conforming
|
to the left initial vector. This vector is overwritten.
|
to the left initial vector. This vector is overwritten.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSGetLeftInitialVector()
|
.seealso: EPSGetLeftInitialVector()
|
|
|
@*/
|
@*/
|
PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,Vec vec)
|
PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,Vec vec)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscTruth breakdown;
|
PetscTruth breakdown;
|
PetscReal norm;
|
PetscReal norm;
|
Vec w;
|
Vec w;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
|
PetscValidHeaderSpecific(vec,VEC_COOKIE,3);
|
PetscValidHeaderSpecific(vec,VEC_COOKIE,3);
|
|
|
/* For the first step, use the initial vector, otherwise a random one */
|
/* For the first step, use the initial vector, otherwise a random one */
|
if (i==0) {
|
if (i==0) {
|
w = eps->vec_initial_left;
|
w = eps->vec_initial_left;
|
}
|
}
|
else {
|
else {
|
ierr = VecDuplicate(eps->vec_initial_left,&w);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->vec_initial_left,&w);CHKERRQ(ierr);
|
ierr = SlepcVecSetRandom(w);CHKERRQ(ierr);
|
ierr = SlepcVecSetRandom(w);CHKERRQ(ierr);
|
}
|
}
|
|
|
/* Force the vector to be in the range of OP */
|
/* Force the vector to be in the range of OP */
|
ierr = STApplyTranspose(eps->OP,w,vec);CHKERRQ(ierr);
|
ierr = STApplyTranspose(eps->OP,w,vec);CHKERRQ(ierr);
|
|
|
/* Orthonormalize the vector with respect to previous vectors */
|
/* Orthonormalize the vector with respect to previous vectors */
|
ierr = IPOrthogonalize(eps->ip,i,PETSC_NULL,eps->W,vec,PETSC_NULL,&norm,&breakdown,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(eps->ip,i,PETSC_NULL,eps->W,vec,PETSC_NULL,&norm,&breakdown,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
if (breakdown) {
|
if (breakdown) {
|
if (i==0) { SETERRQ(1,"Left initial vector is zero"); }
|
if (i==0) { SETERRQ(1,"Left initial vector is zero"); }
|
else { SETERRQ(1,"Unable to generate more left start vectors"); }
|
else { SETERRQ(1,"Unable to generate more left start vectors"); }
|
}
|
}
|
ierr = VecScale(vec,1/norm);CHKERRQ(ierr);
|
ierr = VecScale(vec,1/norm);CHKERRQ(ierr);
|
|
|
if (i!=0) {
|
if (i!=0) {
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
}
|
}
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|