/*
This implements the subspace iteration method for finding the
eigenpairs associtated to the eigenvalues with largest magnitude.
*/
#include "src/eps/epsimpl.h"
typedef struct {
int inner;
} EPS_SUBSPACE;
#undef __FUNCT__
#define __FUNCT__ "EPSSetUp_SUBSPACE"
static int EPSSetUp_SUBSPACE(EPS eps)
{
int ierr;
PetscFunctionBegin;
ierr = EPSDefaultGetWork(eps,1);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSSetDefaults_SUBSPACE"
static int EPSSetDefaults_SUBSPACE(EPS eps)
{
int ierr, N;
EPS_SUBSPACE *subspace = (EPS_SUBSPACE *)eps->data;
PetscFunctionBegin;
ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
if (eps->ncv) {
if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");
}
else eps->ncv = PetscMax(2*eps->nev,eps->nev+8);
if (!eps->max_it) eps->max_it = PetscMax(100,N);
if (!subspace->inner) {
if (eps->ishermitian) subspace->inner = 10;
else subspace->inner = 4;
}
if (!eps->tol) eps->tol = 1.e-7;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSSolve_SUBSPACE"
static int EPSSolve_SUBSPACE(EPS eps,int *its)
{
int ierr, i, j, k, maxit=eps->max_it, ncv = eps->ncv;
Vec w;
PetscReal relerr, tol=eps->tol;
PetscScalar alpha, *H, *S;
EPS_SUBSPACE *subspace = (EPS_SUBSPACE *)eps->data;
PetscFunctionBegin;
w = eps->work[0];
ierr = PetscMalloc(ncv*ncv*sizeof(PetscScalar),&H);CHKERRQ(ierr);
ierr = PetscMalloc(ncv*ncv*sizeof(PetscScalar),&S);CHKERRQ(ierr);
/* Build Z from initial vector */
ierr = VecCopy(eps->vec_initial,eps->V[0]);CHKERRQ(ierr);
for (k=1; k<ncv; k++) {
ierr = STApply(eps->OP,eps->V[k-1],eps->V[k]); CHKERRQ(ierr);
}
/* QR-Factorize V R = Z */
ierr = EPSQRDecomposition(eps,0,ncv,PETSC_NULL,ncv);CHKERRQ(ierr);
i = 0;
*its = 0;
while (*its<maxit) {
/* Y = OP^inner V */
for (k=i;k<ncv;k++) {
for (j=i;j<subspace->inner;j++) {
ierr = STApply(eps->OP,eps->V[k],w);CHKERRQ(ierr);
ierr = VecCopy(w,eps->V[k]);CHKERRQ(ierr);
}
}
/* QR-Factorize V R = Y */
ierr = EPSQRDecomposition(eps,i,ncv,PETSC_NULL,ncv);CHKERRQ(ierr);
/* compute the projected matrix, H = V^* A V */
for (j=i;j<ncv;j++) {
ierr = STApply(eps->OP,eps->V[j],w);CHKERRQ(ierr);
for (k=i;k<ncv;k++) {
ierr = VecDot(w,eps->V[k],H+(k-i)+(ncv-i)*(j-i));CHKERRQ(ierr);
}
}
/* solve the reduced problem, compute the
eigendecomposition H = S Theta S^* */
ierr = EPSDenseNHEP(ncv-i,H,eps->eigr+i,eps->eigi+i,S);CHKERRQ(ierr);
/* update V = V S */
ierr = EPSReverseProjection(eps,i,ncv-i,S);CHKERRQ(ierr);
/* check eigenvalue convergence */
for (j=i;j<ncv;j++) {
ierr = STApply(eps->OP,eps->V[j],w);CHKERRQ(ierr);
alpha = -eps->eigr[j];
ierr = VecAXPY(&alpha,eps->V[j],w);CHKERRQ(ierr);
ierr = VecNorm(w,NORM_2,&relerr);CHKERRQ(ierr);
eps->errest[j] = relerr;
}
/* lock converged Ritz pairs */
eps->nconv = i;
for (j=i;j<ncv;j++) {
if (eps->errest[j]<tol) {
if (j>eps->nconv) {
ierr = EPSSwapEigenpairs(eps,eps->nconv,j);CHKERRQ(ierr);
}
eps->nconv = eps->nconv + 1;
}
}
i = eps->nconv;
*its = *its + 1;
EPSMonitorEstimates(eps,*its,eps->nconv,eps->errest,ncv);
EPSMonitorValues(eps,*its,eps->nconv,eps->eigr,PETSC_NULL,ncv);
if (eps->nconv>=eps->nev) break;
}
ierr = PetscFree(H);CHKERRQ(ierr);
ierr = PetscFree(S);CHKERRQ(ierr);
if( *its==maxit ) *its = *its - 1;
eps->its = *its;
if( eps->nconv == eps->nev ) eps->reason = EPS_CONVERGED_TOL;
else eps->reason = EPS_DIVERGED_ITS;
#if defined(PETSC_USE_COMPLEX)
for (i=0;i<eps->nconv;i++) eps->eigi[i]=0.0;
#endif
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSView_SUBSPACE"
static int EPSView_SUBSPACE(EPS eps,PetscViewer viewer)
{
EPS_SUBSPACE *subspace = (EPS_SUBSPACE *) eps->data;
int ierr;
PetscTruth isascii;
PetscFunctionBegin;
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
if (!isascii) {
SETERRQ1(1,"Viewer type %s not supported for EPSSUBSPACE",((PetscObject)viewer)->type_name);
}
ierr = PetscViewerASCIIPrintf(viewer,"number of inner iterations: %d\n",subspace->inner);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSSetFromOptions_SUBSPACE"
static int EPSSetFromOptions_SUBSPACE(EPS eps)
{
EPS_SUBSPACE *subspace = (EPS_SUBSPACE *)eps->data;
int ierr,val;
PetscTruth flg;
PetscFunctionBegin;
ierr = PetscOptionsHead("SUBSPACE options");CHKERRQ(ierr);
val = subspace->inner;
ierr = PetscOptionsInt("-eps_subspace_inner","Number of inner iterations","EPSSubspaceSetInner",val,&val,&flg);CHKERRQ(ierr);
if (flg) {ierr = EPSSubspaceSetInner(eps,val);CHKERRQ(ierr);}
ierr = PetscOptionsTail();CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "EPSSubspaceSetInner_SUBSPACE"
int EPSSubspaceSetInner_SUBSPACE(EPS eps,int val)
{
EPS_SUBSPACE *subspace;
PetscFunctionBegin;
subspace = (EPS_SUBSPACE *) eps->data;
subspace->inner = val;
PetscFunctionReturn(0);
}
EXTERN_C_END
#undef __FUNCT__
#define __FUNCT__ "EPSSubspaceSetInner"
/*@
EPSSubspaceSetInner - Sets the number of inner iterations performed by
the Subspace Iteration method.
Collective on EPS
Input Parameters:
+ eps - the eigenproblem solver context
- val - number of iterations to perform
Options Database Key:
. -eps_subspace_inner - Sets the value of the inner iterations
Notes:
This value specifies how many matrix-vector products have to be carried
out in the Subspace Iteration method between the orthogonalization steps.
The default value is 10 for Hermitian problems and 4 for non-Hermitian ones.
Level: advanced
@*/
int EPSSubspaceSetInner(EPS eps,int val)
{
int ierr, (*f)(EPS,int);
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_COOKIE);
ierr = PetscObjectQueryFunction((PetscObject)eps,"EPSSubspaceSetInner_C",(void (**)(void))&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(eps,val);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "EPSCreate_SUBSPACE"
int EPSCreate_SUBSPACE(EPS eps)
{
EPS_SUBSPACE *subspace;
int ierr;
PetscFunctionBegin;
ierr = PetscNew(EPS_SUBSPACE,&subspace);CHKERRQ(ierr);
PetscMemzero(subspace,sizeof(EPS_SUBSPACE));
PetscLogObjectMemory(eps,sizeof(EPS_SUBSPACE));
eps->data = (void *) subspace;
eps->ops->setup = EPSSetUp_SUBSPACE;
eps->ops->setdefaults = EPSSetDefaults_SUBSPACE;
eps->ops->solve = EPSSolve_SUBSPACE;
eps->ops->destroy = EPSDefaultDestroy;
eps->ops->setfromoptions = EPSSetFromOptions_SUBSPACE;
eps->ops->view = EPSView_SUBSPACE;
subspace->inner = 0;
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSSubspaceSetInner_C","EPSSubspaceSetInner_SUBSPACE",EPSSubspaceSetInner_SUBSPACE);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_END