Hi all, I am currently developing a code in C++ and I try to use SLEPC in
conjunction with PETSC to solve an eigenvalue problem. Considering I have a set
of similar matrices, my strategy was to solve one problem, and then reuse the
starting vector to accelerate the resolution of the next problems. I tested
this strategy with a rudimentary code to validate the process, but it does not
seems to work. I am afraid I do not fully understand some of the
particularities of SLEPC or PETSC. The main.cpp fille of my code goes like
this: (see at the bottom of this message). You can see that I create a matrix
(A) and two eigensolver contexts (eps and eps2). First a try to solve the
problem using eps. I then recall the invariant subspace and use it as a
starting vector in eps2. This code seems to work, but eps2 takes systematically
more time to converge than eps. Is there something wrong with this code or my
understanding of how SLEPC/PETSC works? Thank you,
Marc
#include "slepceps.h"
int main(int argc, char *argv[])
{
static char help[] = "Simple Hello World example program in SLEPc\n";
Mat A; /* problem matrix */
EPS eps; /* eigenproblem solver context */
EPS eps2;
EPSType type;
PetscReal error,tol,re,im;
PetscScalar kr,ki;
Vec xr,xi;
PetscInt n=30,i,Istart,Iend,nev,maxit,its,nconv;
PetscErrorCode ierr;
PetscInt a,b,c;
SlepcInitialize(&argc,&argv,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
/*
Create matrix
*/
ierr =
MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A);
CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
for (i=Istart;i<Iend;i++) {
if (i>0) { ierr =
MatSetValue(A,i,i-1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if (i<n-1) { ierr =
MatSetValue(A,i,i+1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
ierr = MatSetValue(A,i,i,2.0,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create eigensolver context
*/
ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
its = 4000;
ierr = EPSSetTolerances(eps,PETSC_DEFAULT,its);CHKERRQ(ierr);
nev = (int)n/2;
EPSSetType(eps,EPSKRYLOVSCHUR);
EPSSetDimensions(eps,nev,n,n);
EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE);
EPSSetFromOptions(eps);
ierr = EPSSolve(eps);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
PetscInt nbr = nconv;
Vec subspace;
VecCreate(PETSC_COMM_WORLD,&subspace);
VecSetSizes(subspace,PETSC_DECIDE,n);
VecSetFromOptions(subspace);
Vec *sub;
VecDuplicateVecs(subspace,nbr,&sub);
EPSGetInvariantSubspace(eps,sub);
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = EPSCreate(PETSC_COMM_WORLD,&eps2);CHKERRQ(ierr);
ierr = EPSSetInitialSpace(eps2,nbr,sub);
ierr = EPSSetOperators(eps2,A,NULL);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps2,EPS_HEP);CHKERRQ(ierr);
its = 4000;
ierr = EPSSetTolerances(eps2,PETSC_DEFAULT,its);CHKERRQ(ierr);
nev = (int)n/2;
ierr = EPSSetType(eps2,EPSKRYLOVSCHUR);CHKERRQ(ierr);
ierr = EPSSetDimensions(eps2,nev,n,n);CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps2,EPS_SMALLEST_MAGNITUDE);CHKERRQ(ierr);
ierr = EPSSetFromOptions(eps2);CHKERRQ(ierr);
ierr = EPSSolve(eps2);CHKERRQ(ierr);
// ierr = EPSGetConverged(eps2,&nconv);CHKERRQ(ierr);
ierr = EPSDestroy(&eps2);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);
ierr = VecDestroy(&xi);CHKERRQ(ierr);
ierr = SlepcFinalize();
return 0;
}
