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

[Matthew Knepley]

On Wed, May 28, 2014 at 12:27 PM, Giacomo Mulas
<gmulas at oa-cagliari.inaf.it>wrote:

> 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?
>
> 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?
>
> 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?
>
> 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?
>

There is not. However, since this is just a set of Vec pointers, allocating
and copying
should be fine. The amount of memory taken up by the pointers is very very
small.

  Thanks,

     Matt


> 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)
> _________________________________________________________________
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140528/68a50eb5/attachment.html>