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