<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, May 28, 2014 at 12:27 PM, Giacomo Mulas <span dir="ltr"><<a href="mailto:gmulas@oa-cagliari.inaf.it" target="_blank">gmulas@oa-cagliari.inaf.it</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello.<br>
<br>
After some time stuck doing less exciting stuff, I got back to trying to use<br>
slepc for science. I am trying to use the relatively new functionality for<br>
arbitrary selection of the eigenpairs to be found, and I would want some<br>
support to understand if I am doing things correctly. I apologise with Jose<br>
Roman for disappearing when I should have helped testing this, it was not my<br>
choice. But now I am back at it.<br>
<br>
In particular, I want to obtain the eigenvectors with the maximum projection<br>
in a given subspace (defined by a number of normalised vectors, let's call<br>
them targets).<br>
<br>
I don't know in advance how many eigenpairs must be determined: I want to<br>
obtain enough eigenvectors that the projection of all targets in the space<br>
of these eigenvectors is very nearly identical to the targets themselves.<br>
<br>
So my strategy, so far, is the following:<br>
<br>
1) create and set up the matrix H to be diagonalised<br>
<br>
ierr = MatCreate(mixpars->slepc_comm, &H); CHKERRQ(ierr);<br>
ierr = MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, statesinlist,<br>
statesinlist);<br>
CHKERRQ(ierr);<br>
ierr = MatSetFromOptions(H);CHKERRQ(<u></u>ierr);<br>
ierr = MatSetUp(H);CHKERRQ(ierr);<br>
...<br>
... some MatSetValue(H,...);<br>
...<br>
ierr = MatAssemblyBegin(H,MAT_FINAL_<u></u>ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(H,MAT_FINAL_<u></u>ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatGetVecs(H, &xr, PETSC_NULL);CHKERRQ(ierr);<br>
ierr = MatGetVecs(H, &xi, PETSC_NULL);CHKERRQ(ierr);<br>
<br>
2) create the eps<br>
<br>
ierr = EPSCreate(mixpars->slepc_comm,<u></u>&eps);CHKERRQ(ierr);<br>
ierr = EPSSetOperators(eps,H,PETSC_<u></u>NULL);CHKERRQ(ierr);<br>
ierr = EPSSetProblemType(eps,EPS_HEP)<u></u>;CHKERRQ(ierr);<br>
ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE);<br>
ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);<br>
ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);<br>
ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);<br>
ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);<br>
<br>
3) set the eps to use an arbitrary selection function<br>
<br>
/* every time I solve, I want to find one eigenvector, and it must be<br>
the one with the largest component along the target state */<br>
ierr = EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE);<br>
ierr = EPSSetArbitrarySelection(eps, computeprojection(),<br>
(void *) &targetcompindex);<br>
ierr = EPSSetDimensions(eps, 1, PETSC_IGNORE, PETSC_IGNORE);<br>
CHKERRQ(ierr);<br>
<br>
4) run a loop over the target vectors, and iteratively call EPSSolve until<br>
each target vector is completely contained in the space of the<br>
eigenvectors found. Before each call to EPSSolve, I set the initial guess<br>
equal to the target vector, and set the deflation space to be the set of<br>
eigenvectors found so far. After each call to EPSSolve, I add the new<br>
eigenvectors to the deflation space one by one, and check if the target<br>
state is (nearly) fully contained in the eigenvectors space. If yes, I<br>
move on to the next target state and so on.<br>
<br>
5) free everything, destroy eps, matrices, vectors etc.<br>
<br>
<br>
I have some questions about the above:<br>
<br>
1) should it work in principle, or am I getting it all wrong?<br>
<br>
2) should I destroy and recreate the eps after each call to EPSSolve and<br>
before next call? Or, since the underlying matrix is always the same, can I<br>
just call EPSSetInitialSpace(), EPSSetDeflationSpace(), update the internal<br>
parameter to be passed to the arbitrary selection function and I can call<br>
again EPSSolve?<br>
<br>
3) Since what I want is going on to find one eigenpair at a time of the same<br>
problem until some condition is fulfilled, is there a way in which I can<br>
achieve this without setting it up again and again every time? Can I specify<br>
an arbitrary function that is called by EPSSolve to decide whether enough<br>
eigenpairs were computed or not, instead of doing it in this somewhat<br>
awkward manner?<br>
<br>
4) more technical: since I add vectors one by one to the deflation space, to<br>
begin with I allocate a Vec *Cv with PetscMalloc(statesinlist*<u></u>sizeof(Vec *), &Cv);<br>
where statesinlist is the size of the problem, hence the maximum<br>
hypothetical size of the deflation space. I would prefer to allocate this<br>
dynamically, enlarging it as needed. Is there something like realloc() in<br>
PETSC/SLEPC?<br></blockquote><div><br></div><div>There is not. However, since this is just a set of Vec pointers, allocating and copying</div><div>should be fine. The amount of memory taken up by the pointers is very very small.</div>
<div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Thanks in advance, bye<br>
Giacomo Mulas<br>
<br>
-- <br>
______________________________<u></u>______________________________<u></u>_____<br>
<br>
Giacomo Mulas <<a href="mailto:gmulas@oa-cagliari.inaf.it" target="_blank">gmulas@oa-cagliari.inaf.it</a>><br>
______________________________<u></u>______________________________<u></u>_____<br>
<br>
INAF - Osservatorio Astronomico di Cagliari<br>
via della scienza 5 - 09047 Selargius (CA)<br>
<br>
tel. <a href="tel:%2B39%20070%2071180244" value="+3907071180244" target="_blank">+39 070 71180244</a><br>
mob. : <a href="tel:%2B39%20329%20%206603810" value="+393296603810" target="_blank">+39 329 6603810</a><br>
______________________________<u></u>______________________________<u></u>_____<br>
<br>
"When the storms are raging around you, stay right where you are"<br>
(Freddy Mercury)<br>
______________________________<u></u>______________________________<u></u>_____<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>