[petsc-users] question about arbitrary eigenvector selection in SLEPC

[Jose E. Roman]

El 28/05/2014, a las 19:27, Giacomo Mulas escribió:

> Hello.
> 
> After some time stuck doing less exciting stuff, I got back to trying to use
> slepc for science.  I am trying to use the relatively new functionality for
> arbitrary selection of the eigenpairs to be found, and I would want some
> support to understand if I am doing things correctly.  I apologise with Jose
> Roman for disappearing when I should have helped testing this, it was not my
> choice. But now I am back at it.
> 
> In particular, I want to obtain the eigenvectors with the maximum projection
> in a given subspace (defined by a number of normalised vectors, let's call
> them targets).
> 
> I don't know in advance how many eigenpairs must be determined: I want to
> obtain enough eigenvectors that the projection of all targets in the space
> of these eigenvectors is very nearly identical to the targets themselves.
> 
> So my strategy, so far, is the following:
> 
> 1) create and set up the matrix H to be diagonalised
> 
>  ierr = MatCreate(mixpars->slepc_comm, &H); CHKERRQ(ierr);
>  ierr = MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, statesinlist,
>                     statesinlist);
>  CHKERRQ(ierr);
>  ierr = MatSetFromOptions(H);CHKERRQ(ierr);
>  ierr = MatSetUp(H);CHKERRQ(ierr);
>  ...
>  ... some MatSetValue(H,...);
>  ...
>  ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>  ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>  ierr = MatGetVecs(H, &xr, PETSC_NULL);CHKERRQ(ierr);
>  ierr = MatGetVecs(H, &xi, PETSC_NULL);CHKERRQ(ierr);
> 
> 2) create the eps
> 
>  ierr = EPSCreate(mixpars->slepc_comm,&eps);CHKERRQ(ierr);
>  ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr);
>  ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
>  ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE);
>  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
>  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
>  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
>  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
> 
> 3) set the eps to use an arbitrary selection function
> 
>  /* every time I solve, I want to find one eigenvector, and it must be
>     the one with the largest component along the target state */
>  ierr = EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE);
>  ierr = EPSSetArbitrarySelection(eps, computeprojection(),
> 				  (void *) &targetcompindex);
>  ierr = EPSSetDimensions(eps, 1, PETSC_IGNORE, PETSC_IGNORE);
>  CHKERRQ(ierr);
> 
> 4) run a loop over the target vectors, and iteratively call EPSSolve until
>   each target vector is completely contained in the space of the
>   eigenvectors found. Before each call to EPSSolve, I set the initial guess
>   equal to the target vector, and set the deflation space to be the set of
>   eigenvectors found so far. After each call to EPSSolve, I add the new
>   eigenvectors to the deflation space one by one, and check if the target
>   state is (nearly) fully contained in the eigenvectors space. If yes, I
>   move on to the next target state and so on.
> 
> 5) free everything, destroy eps, matrices, vectors etc.
> 
> 
> I have some questions about the above:
> 
> 1) should it work in principle, or am I getting it all wrong?

I don't see much problem.

> 
> 2) should I destroy and recreate the eps after each call to EPSSolve and
> before next call? Or, since the underlying matrix is always the same, can I
> just call EPSSetInitialSpace(), EPSSetDeflationSpace(), update the internal
> parameter to be passed to the arbitrary selection function and I can call
> again EPSSolve?

No need to recreate the solver. The only thing is EPSSetDeflationSpace() - I would suggest calling EPSRemoveDeflationSpace() and then EPSSetDeflationSpace() again with the extended set of vectors. Do not call EPSSetDeflationSpace() with a single vector every time.

> 
> 3) Since what I want is going on to find one eigenpair at a time of the same
> problem until some condition is fulfilled, is there a way in which I can
> achieve this without setting it up again and again every time? Can I specify
> an arbitrary function that is called by EPSSolve to decide whether enough
> eigenpairs were computed or not, instead of doing it in this somewhat
> awkward manner?

No.

> 
> 4) more technical: since I add vectors one by one to the deflation space, to
> begin with I allocate a Vec *Cv with PetscMalloc(statesinlist*sizeof(Vec *), &Cv);
> where statesinlist is the size of the problem, hence the maximum
> hypothetical size of the deflation space. I would prefer to allocate this
> dynamically, enlarging it as needed. Is there something like realloc() in
> PETSC/SLEPC?
> 
> Thanks in advance, bye
> Giacomo Mulas
> 
> -- 
> _________________________________________________________________
> 
> Giacomo Mulas <gmulas at oa-cagliari.inaf.it>
> _________________________________________________________________
> 
> INAF - Osservatorio Astronomico di Cagliari
> via della scienza 5 - 09047 Selargius (CA)
> 
> tel.   +39 070 71180244
> mob. : +39 329  6603810
> _________________________________________________________________
> 
> "When the storms are raging around you, stay right where you are"
>                         (Freddy Mercury)
> _________________________________________________________________