[petsc-users] question about arbitrary eigenvector selection in SLEPC
Jose E. Roman
jroman at dsic.upv.es
Wed May 28 13:59:27 CDT 2014
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)
> _________________________________________________________________
More information about the petsc-users
mailing list